Explaining the Variability in the Profiles of Ratings of Perceived Exertion for a Dynamic Upper Extremity Task: A Functional Regression Approach

Authors

Affiliations

Setareh Kazemi Kheiri

Department of Industrial and Systems Engineering, University at Buffalo

Sahand Hajifar

Department of Industrial and Systems Engineering, University at Buffalo

Linh Tran

Farmer School of Business, Miami University

Zahra Vahedi

Department of Industrial and Systems Engineering, University at Buffalo

Lora A. Cavuoto

Department of Industrial and Systems Engineering, University at Buffalo

Fadel M. Megahed

Farmer School of Business, Miami University

Hongyue Sun

College of Engineering, University of Georgia

Objectives of this Document

In this document all the code, results and analysis related to the associated paper entitled “Explaining the Variability in the Profiles of Ratings of Perceived Exertion for a Dynamic Upper Extremity Task: A Functional Regression Approach” are provided. The main objective of doing this research was gaining a better understanding of how different factors contribute to the development of fatigue during manual material handling (MMH) operations in a simulated warehousing environment. Our approach is divided into the following main steps:

1. Data pre-processing;

2. Data transformation;

3. Functional feature extraction;

4. Functional regression on the transformed data with different set of scalar and functional features;

5. Benchmark models;

6. Checking models’ assumptions (residual analysis)

---

1 R Setup and Required Packages

In this project, the open-source programming language is used for our analysis. is maintained by an international team of developers who make the language available at The Comprehensive R Archive Network. Readers interested in reusing our code and reproducing our results should have installed locally on their machines. can be installed on a number of different operating systems (see Windows, Mac, and Linux for the installation instructions for these systems). We also recommend using the RStudio interface for . The reader can download RStudio for free by following the instructions at the link. For non-R users, we recommend the Hands-on Programming with R for a brief overview of the software’s functionality. Hereafter, we assume that the reader has an introductory understanding of the programming language.

# Setting the random seed and chunk dependencies
knitr::opts_chunk$set(cache.extra = set.seed(2022),
                      autodep = TRUE) 

# Load required packages
if(require(pacman)==F) install.packages('pacman')
pacman::p_load(tidyverse, stringr, magrittr, purrr, imputeTS, timetk, refund, forecast, zoo, DT, stargazer, readxl, grid, gridExtra, fda)

# Graphic setup
theme_set(theme_bw() +
          theme(legend.position = 'top',
                axis.title = element_text(size=8)))
scale_colour_discrete = scale_color_brewer(palette = "Dark2")

# Load features data 
load("Features.RData")

# Load TRPE 
load("y_T_val.RData")

# Load RPE
load("y.RData")

# Load Anthropometric data
load("Anthrpmtrc_data.RData")

# Load raw IMU data
# load("IMU_raw.RData")

# Load the processed data from IMU_raw.RData due to the large file, solely for the purpose of rendering the QMD file
load("sensors.RData")

design = readxl::read_excel("Experiment Design.xlsx", col_names = F) #The excel file which shows the random order of task conditions assigned to each subject in each session of the experiment

CP = read.csv("corrected_Changepoints.csv") #The changepoints based on the raw sensors input

2 Data Preprocessing

2.1 Overview

2.1.1 Features

# We excluded FFT features for the purpose of this study
features = features_all |> 
  dplyr::select(-c(dplyr::contains("FFT"), cycle, Acc_Length))

# Display features
DT::datatable(
  data = features, rownames = F, 
  extensions = c('Buttons','FixedColumns'),
  options = list(
    dom = 'Bfrtip',
    buttons = c('csv', 'excel', 'pdf', 'print'),
    paging = TRUE, searching = TRUE, info = FALSE,
    sort = TRUE, scrollX = TRUE, fixedColumns = list(leftColumns = 3)
  )
) |> 
  DT::formatRound(
    columns = 5:ncol(features), 
    digits = 3)

2.1.2 TRPE

The TRPE is the transformed RPE, using a multivariate Box-Cox transformation. Details on how this was implemented can be found in our previous work.

In the code chunk below, we assigned the row names and column names for the TRPE data (i.e. y_T_val) as those from the orignial RPE data (i.e. y) since the two shared the same structure.

## extract subject-condition from the RPE df and reformat
## e.g. from "Sub01 Condition 1.5-15" to "1-1.5-15"
y_sub_con = rownames(y) |> 
  stringr::str_c(collapse = ",") |> 
  stringr::str_replace_all(c(" Condition " = "-", "Sub(0)?" = "")) |> 
  stringr::str_split(pattern = ",") |> 
  magrittr::extract2(1)

## assign col names and row names for y_T_val
rownames(y_T_val) = y_sub_con
colnames(y_T_val) = as.numeric(substr(colnames(y),2,3)) 

## Display TRPE
DT::datatable(
  data = y_T_val, 
  extensions = c('Buttons','FixedColumns'),
  options = list(
    dom = 'Bfrtip',
    buttons = c('csv', 'excel', 'pdf', 'print'),
    paging = TRUE, searching = TRUE, info = FALSE,
    sort = TRUE, scrollX = TRUE, fixedColumns = list(leftColumns = 1)
  )
) |> 
  DT::formatRound(
    columns = 1:ncol(y_T_val), 
    digits = 3)

2.2 Prepare the input data for modeling

For the purpose of this study, we only selected features of those cycles whose median time was the closest to the time when an RPE was recorded.

It must also be noted that the extracted features contained 50 subject-sessions after removing those whose features were considered outliers, while the TRPE data contained 44, wherein the participants successfully completed the full first 45 minutes. As a result, it is important to only select the data of the subject-sessions whose features and TRPE are both available as the input for our models.

features = features |> 
  # filter out NAs (i.e. only select the desired cycles)
  dplyr::filter(!is.na(rpe)) |> 
  dplyr::mutate(sub_con = str_c(subject, str_split(condition, " ", simplify = T)[,2], sep = "-"), .keep="unused", .before = Time,
                subject = as.character(subject)) |> 
  # only keep data of the matched subject-sessions (or subject-conditions)
  dplyr::filter(sub_con %in% y_sub_con) |> 
  dplyr::group_by(sub_con) |> 
  dplyr::arrange(Time, .by_group = T) |> 
  dplyr::mutate(Time = as.integer(5*(row_number() - 1))) |> 
  dplyr::select(-c(session, rpe)) |> 
  dplyr::ungroup()

In the code chunk below, we can see that there were a total of 37 subject-sessions, each had 10 observations (corresponding to the 10 RPEs recorded every 5 minutes during the 45-minute interval of the experiment). However, subject 10 - session 4 and subject 17 - session 3 only contained 9 observations. This was because the cycle associated with their 10th RPE was considered an outlier and was removed during the feature extraction step.

summary = features |> 
  dplyr::group_by(sub_con) |> 
  dplyr::tally()

DT::datatable(
  data = summary, rownames = F, 
  extensions = c('Buttons','FixedColumns'),
  options = list(
    dom = 'Bfrtip',
    buttons = c('csv', 'excel', 'pdf', 'print'),
    paging = TRUE, searching = TRUE, info = FALSE,
    sort = TRUE, scrollX = TRUE)
  ) 

To overcome this issue, given that there were only two subject-sessions that missed the features for one time point (out of 10), we decided to impute the missing values to make use of the data that we got.

In the code chunk below, we compare several imputation methods (median, naive, simple moving average, weighted moving average, and Kalman smoothing) to see which one results in the lowest average MAE, MAPE, RMSE across all features, using the subject-sessions that have full data.

# get subject-sessions with full 10 RPEs

full_rpe = features |> 
  dplyr::filter(!sub_con %in% c("10-2.5-5", "17-2.5-5")) |> 
  dplyr::select(-c(subject))

# try different imputation methods
median = full_rpe |> 
  dplyr::group_by(sub_con) |> 
  dplyr::summarise(across(-c("Time"), \(x) median(x))) |> 
  dplyr::select(-c(sub_con))

naive = full_rpe |> 
  dplyr::group_by(sub_con) |> 
  dplyr::filter(Time == 40) |> 
  dplyr::ungroup() |> 
  dplyr::select(-c(Time, sub_con))

sma = full_rpe |> 
  dplyr::group_by(sub_con) |> 
  dplyr::mutate(across(-c("Time"), \(x) if_else(Time == 45, NA_real_, x))) |> 
  dplyr::mutate(across(-c("Time"), \(x) imputeTS::na_ma(x, k=3, weighting = "simple"))) |> 
  dplyr::filter(Time == 45) |> 
  dplyr::ungroup() |>
  dplyr::select(-c(Time, sub_con))

lwma = full_rpe |> 
  dplyr::group_by(sub_con) |> 
  dplyr::mutate(across(-c("Time"), \(x) if_else(Time == 45, NA_real_, x))) |> 
  dplyr::mutate(across(-c("Time"), \(x) imputeTS::na_ma(x, k=3, weighting = "linear"))) |> 
  dplyr::filter(Time == 45) |> 
  dplyr::ungroup() |>
  dplyr::select(-c(Time, sub_con))

ewma = full_rpe |> 
  dplyr::group_by(sub_con) |> 
  dplyr::mutate(across(-c("Time"), \(x) if_else(Time == 45, NA_real_, x))) |> 
  dplyr::mutate(across(-c("Time"), \(x) imputeTS::na_ma(x, k=3))) |> 
  dplyr::filter(Time == 45) |> 
  dplyr::ungroup() |>
  dplyr::select(-c(Time, sub_con))
  

kalman = full_rpe |> 
  dplyr::group_by(sub_con) |> 
  dplyr::mutate(across(-c("Time"), \(x) if_else(Time == 45, NA_real_, x))) |> 
  dplyr::mutate(across(-c("Time"), \(x) na_kalman(x))) |> 
  dplyr::filter(Time == 45) |> 
  dplyr::ungroup() |>
  dplyr::select(-c(Time, sub_con))

# compare

org = full_rpe |> 
  dplyr::filter(Time == 45) |> 
  dplyr::select(-c(sub_con, Time))

metrics_extract = function(imp, org) {
  metrics = purrr::map2_dfr(.x =imp, .y = org, .f = forecast::accuracy)

MAE = metrics |>
  dplyr::mutate(across(1:54, \(x) magrittr::extract2(x, 3))) |> 
  tidyr::pivot_longer(everything(), names_to = "features", values_to = "MAE")

MAPE = metrics |>
  dplyr::mutate(across(1:54, \(x) magrittr::extract2(x, 5))) |> 
  tidyr::pivot_longer(everything(), names_to = "features", values_to = "MAPE")

RMSE = metrics |>
  dplyr::mutate(across(1:54, \(x) magrittr::extract2(x, 2))) |> 
  tidyr::pivot_longer(everything(), names_to = "features", values_to = "RMSE")

all_metrics = MAE |> 
  dplyr::inner_join(MAPE, by = "features") |> 
  dplyr::inner_join(RMSE, by = "features") |> 
  dplyr::summarise(across(c("MAE", "MAPE", "RMSE"),
                   list(Avg = mean)))
}


summary = dplyr::bind_rows(
  metrics_extract(median, org) |> dplyr::mutate(method = "Median"),
  metrics_extract(naive, org) |> dplyr::mutate(method = "Naive"),
  metrics_extract(sma, org) |> dplyr::mutate(method = "SMA"),
  metrics_extract(lwma, org) |> dplyr::mutate(method = "LWMA"),
  metrics_extract(ewma, org) |> dplyr::mutate(method = "EWMA"),
  metrics_extract(kalman, org) |> dplyr::mutate(method = "Kalman")) |> 
  dplyr::relocate(method, .before = "MAE_Avg")

colnames(summary) = c('Methods', 'Avg MAE', 'Avg MAPE', 'Avg RMSE')

DT::datatable(
  data = summary, rownames = F, 
  options = list(
    dom = 'Bfrtip',
    info = FALSE,
    sort = TRUE, scrollX = TRUE)
  ) 

2.2.1 Imputation

Based on the result above, median is used to impute the missing features.

# Add 2 rows of NA feature values for subcon "10-2.5-5" and "17-2.5-5" at time 45
na_rows = dplyr::bind_rows(tibble::as_tibble_row(c(10, "10-2.5-5", 45, rep(NA, 54)), .name_repair = 'unique'),
                    tibble::as_tibble_row(c(17, "17-2.5-5", 45, rep(NA, 54)), .name_repair = 'unique'))
colnames(na_rows) = colnames(features)

na_rows = na_rows |> 
  dplyr::mutate(across(.cols = 4:57, .fns = as.double),
         Time = as.integer(Time))

features = dplyr::bind_rows(features, na_rows)

## Perform imputation with median
features_imputed = features |> 
  dplyr::group_by(sub_con) |> 
  dplyr::mutate(across(.cols = everything(), \(x) replace_na(x, median(x, na.rm=T))))

2.2.2 Joining the data together

We then join the feature data with the TRPE and anthropometric data to get a complete dataset.

Anthrpmtrc_data_c = Anthrpmtrc_data |> 
  dplyr::mutate(subject = stringr::str_extract(`Sub#`, "[0-9]{2}") |> as.integer() |> as.character(), .keep = "unused") |> 
  dplyr::select(-c(Session))

input_data_nested = 
  tibble::as_tibble(y_T_val, rownames = "sub_con") |> 
  tidyr::pivot_longer(cols = -c(1), names_to = "Time", values_to = "TRPE") |> 
  dplyr::mutate(xpace = str_split(sub_con, "-", simplify = T)[,2] |> as.numeric(),
                xload = str_split(sub_con, "-", simplify = T)[,3]|> as.numeric(),
                Time = as.integer(Time)) |> 
  dplyr::inner_join(features_imputed, by = c("sub_con", "Time")) |> 
  dplyr::group_by(subject, sub_con, xpace, xload) |> 
  tidyr::nest() |> 
  dplyr::left_join(Anthrpmtrc_data_c, by = "subject") |> 
  dplyr::mutate(Gender = as.factor(Gender))

colnames(input_data_nested)[6:11] = c("xgender", "xage", "xheight", "xweight", "xwaist", "xhip")

## Extract list-columns

for (i in 1:56) {
  input_data_nested = input_data_nested |> 
    dplyr::mutate(purrr::map(.x = data, .f = magrittr::extract2, i))
  colnames(input_data_nested)[i+11] = colnames(input_data_nested$data[[1]])[[i]]
}

input_data_nested = input_data_nested |> 
  dplyr::ungroup() |> 
  dplyr::select(-c(data, subject))

2.2.3 For Functional Regression Models

Due to the unique input format of the pffr(), which is the function that we used to fit our functional regression models, this section is dedicated to transforming the input data into the suitable format.

subcon = input_data_nested$sub_con ## extract sub-con
time = input_data_nested$Time[[1]] ## extract time (i.e., 0, 5, 10, ... 45)

## Copy scalar variables to an object named `input_data`
input_data = c(input_data_nested[, 2:9])

## Repeat the procedure above to the all features
for (i in 11:65) {
  ## create obj where the column is extracted and transformed to matrix form
  ## make obj a list to append to `input_data`
  obj = matrix(unlist(input_data_nested[i]), nrow = 37, ncol=10, byrow = T, 
               dimnames = list(subcon, time)) |> list()
  
  input_data = c(input_data, obj) ## append
  
  names(input_data)[[i-2]] = colnames(input_data_nested[i]) ## assign name of the element with the corresponding feature name 
}

names(input_data)[[9]] = "Y"

3 Functional Regression

After preparing the inputs, they are fed into the pffr function of refund package. It is important to specify the type of covariates of the model based on the multivariate model specified in their package # multivariate model:E(Y(t)) = _0(t) + X1(s)_1(s,t)ds + xlin _3(t) + f_1(xte1, xte2) + f_2(xsmoo, t) + _4 xconst. Since anthropometric data and task specifications (load and pace) are constant through the experiments, they are considered as constant variables (xconst). The features from sensors are smooth functions varying over y index (time), so they are considered as xsmooth.

3.1 With only task specifications

# Model 1 only considers task specifications
model_1 <- pffr(Y ~
                 c(xload)  +  # linear effect of xload constant over Y-index
                 c(xpace),
                 data=input_data, method = "REML", algorithm = "gam")
summary(model_1)

Family: gaussian 
Link function: identity 

Formula:
Y ~ c(xload) + c(xpace)

Constant coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.68991    0.45888  -1.503    0.134    
xload        0.08898    0.01598   5.570 4.99e-08 ***
xpace        0.63447    0.14873   4.266 2.55e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Smooth terms & functional coefficients:
                    edf Ref.df     F p-value    
Intercept(yindex) 6.228      9 28.87  <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.441   Deviance explained = 45.4%
-REML score = 562.72  Scale est. = 1.1486    n = 370 (37 x 10)

3.2 Adding individual characteristics

# Model 2 considers the individual characteristics in addition to the task specifications
model_2 <- pffr(Y ~
                 c(xload)  +
                 c(xpace)  +  
                 c(xgender)+ 
                 c(xage)+
                 c(xheight)+ 
                 c(xweight)+ 
                 c(xhip)+ 
                 c(xwaist),
                 data=input_data, method = "REML", algorithm = "gam")
summary(model_2)

Family: gaussian 
Link function: identity 

Formula:
Y ~ c(xload) + c(xpace) + c(xgender) + c(xage) + c(xheight) + 
    c(xweight) + c(xhip) + c(xwaist)

Constant coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.662385   1.939223   0.342 0.732875    
xload        0.101437   0.012607   8.046 1.30e-14 ***
xpace        0.789816   0.115773   6.822 3.90e-11 ***
xgenderF    -1.374987   0.412211  -3.336 0.000941 ***
xage         0.086763   0.019867   4.367 1.65e-05 ***
xheight     -0.044392   0.011958  -3.712 0.000238 ***
xweight      0.014993   0.006086   2.463 0.014234 *  
xhip         0.082202   0.009269   8.868  < 2e-16 ***
xwaist      -0.062869   0.011987  -5.245 2.70e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Smooth terms & functional coefficients:
                    edf Ref.df     F p-value    
Intercept(yindex) 6.949      9 48.84  <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.666   Deviance explained =   68%
-REML score = 488.24  Scale est. = 0.68574   n = 370 (37 x 10)

3.3 With only the individual characteristics

# Model 3 considers only the individual characteristics
model_3 <- pffr(Y ~
                 c(xgender)+ 
                 c(xage)+
                 c(xheight)+ 
                 c(xweight)+ 
                 c(xhip)+ 
                 c(xwaist),
                 data=input_data, method = "REML", algorithm = "gam")
summary(model_3)

Family: gaussian 
Link function: identity 

Formula:
Y ~ c(xgender) + c(xage) + c(xheight) + c(xweight) + c(xhip) + 
    c(xwaist)

Constant coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.329940   2.110020   0.630  0.52890    
xgenderF    -0.902932   0.447459  -2.018  0.04435 *  
xage         0.065794   0.021594   3.047  0.00248 ** 
xheight     -0.029231   0.012937  -2.260  0.02445 *  
xweight      0.010712   0.006639   1.613  0.10754    
xhip         0.076002   0.010108   7.519 4.52e-13 ***
xwaist      -0.051474   0.013048  -3.945 9.62e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Smooth terms & functional coefficients:
                    edf Ref.df     F p-value    
Intercept(yindex) 6.703      9 40.53  <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.599   Deviance explained = 61.3%
-REML score = 517.04  Scale est. = 0.82366   n = 370 (37 x 10)

3.4 Adding statistical features from the sensors in addition to the individual characteristics and task specifications

# Below models consider features from sensors in addition to the individual characteristics and task specifications.

model_med <- pffr(Y ~
                 c(xload)  +
                 c(xpace)  +  
                 c(xgender)+ 
                 c(xage)+
                 c(xheight)+ 
                 c(xweight)+ 
                 c(xhip)+ 
                 c(xwaist)+
                 +s( Wrist_Acc_Med )
                 +s( Wrist_Jerk_Med )
                 +s( Trunk_Acc_Med )
                 +s( Trunk_Jerk_Med )
                 +s( Upperarm_Acc_Med )
                 +s( Upperarm_Jerk_Med ),
                 data=input_data, method = "REML", algorithm = "gam")
summary(model_med)

Family: gaussian 
Link function: identity 

Formula:
Y ~ c(xload) + c(xpace) + c(xgender) + c(xage) + c(xheight) + 
    c(xweight) + c(xhip) + c(xwaist) + +s(Wrist_Acc_Med) + s(Wrist_Jerk_Med) + 
    s(Trunk_Acc_Med) + s(Trunk_Jerk_Med) + s(Upperarm_Acc_Med) + 
    s(Upperarm_Jerk_Med)

Constant coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.372477   1.817378   0.205 0.837743    
xload        0.093877   0.026106   3.596 0.000376 ***
xpace        0.873847   0.103965   8.405 1.55e-15 ***
xgenderF    -2.290008   0.381916  -5.996 5.60e-09 ***
xage         0.067665   0.018404   3.677 0.000278 ***
xheight     -0.055371   0.011437  -4.841 2.03e-06 ***
xweight      0.012528   0.006452   1.942 0.053058 .  
xhip         0.091989   0.009127  10.079  < 2e-16 ***
xwaist      -0.038017   0.011455  -3.319 0.001011 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Smooth terms & functional coefficients:
                        edf Ref.df      F p-value    
Intercept(yindex)     7.764  9.000 92.123  <2e-16 ***
s(Wrist_Acc_Med)      9.816 13.775  2.142  0.0104 *  
s(Wrist_Jerk_Med)     2.103  2.905  3.567  0.0176 *  
s(Trunk_Acc_Med)     12.572 16.861  6.168  <2e-16 ***
s(Trunk_Jerk_Med)     6.046  6.671  7.768  <2e-16 ***
s(Upperarm_Acc_Med)   6.153  6.730 11.333  <2e-16 ***
s(Upperarm_Jerk_Med)  5.225  6.361  9.649  <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.827   Deviance explained = 85.4%
-REML score = 409.61  Scale est. = 0.35536   n = 370 (37 x 10)

model_mean <- pffr(Y ~
                 c(xload)  + 
                 c(xpace)  +  
                 c(xgender)+ 
                 c(xage)+
                 c(xheight)+ 
                 c(xweight)+ 
                 c(xhip)+ 
                 c(xwaist)+
                 +s( Wrist_Acc_Mean )
                 +s( Wrist_Jerk_Mean )
                 +s( Trunk_Acc_Mean )
                 +s( Trunk_Jerk_Mean )
                 +s( Upperarm_Acc_Mean )
                 +s( Upperarm_Jerk_Mean ),
                 data=input_data, method = "REML", algorithm = "gam")
summary(model_mean)

Family: gaussian 
Link function: identity 

Formula:
Y ~ c(xload) + c(xpace) + c(xgender) + c(xage) + c(xheight) + 
    c(xweight) + c(xhip) + c(xwaist) + +s(Wrist_Acc_Mean) + s(Wrist_Jerk_Mean) + 
    s(Trunk_Acc_Mean) + s(Trunk_Jerk_Mean) + s(Upperarm_Acc_Mean) + 
    s(Upperarm_Jerk_Mean)

Constant coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.915319   2.055391   0.932 0.352215    
xload        0.075870   0.021467   3.534 0.000478 ***
xpace        0.776461   0.100445   7.730 1.93e-13 ***
xgenderF    -2.080533   0.418962  -4.966 1.19e-06 ***
xage         0.105941   0.020447   5.181 4.22e-07 ***
xheight     -0.070594   0.013171  -5.360 1.74e-07 ***
xweight      0.034454   0.007228   4.767 3.01e-06 ***
xhip         0.096624   0.009672   9.990  < 2e-16 ***
xwaist      -0.057875   0.012184  -4.750 3.25e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Smooth terms & functional coefficients:
                         edf Ref.df      F  p-value    
Intercept(yindex)      7.901  9.000 97.480  < 2e-16 ***
s(Wrist_Acc_Mean)     11.461 15.241  3.100 9.47e-05 ***
s(Wrist_Jerk_Mean)    13.861 18.161  2.656 0.000329 ***
s(Trunk_Acc_Mean)     11.391 15.157  5.499  < 2e-16 ***
s(Trunk_Jerk_Mean)    14.839 19.535  2.551 0.000423 ***
s(Upperarm_Acc_Mean)  17.198 20.726  3.981  < 2e-16 ***
s(Upperarm_Jerk_Mean)  4.021  4.415  3.639 0.003233 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.851   Deviance explained = 88.7%
-REML score = 416.16  Scale est. = 0.30654   n = 370 (37 x 10)

model_Perc25 <- pffr(Y ~
                 c(xload)  +
                 c(xpace)  +  
                 c(xgender)+ 
                 c(xage)+
                 c(xheight)+ 
                 c(xweight)+ 
                 c(xhip)+ 
                 c(xwaist)+
                 +s( Wrist_Acc_Perc25 )
                 +s( Wrist_Jerk_Perc25 )
                 +s( Trunk_Acc_Perc25 )
                 +s( Trunk_Jerk_Perc25 )
                 +s( Upperarm_Acc_Perc25 )
                 +s( Upperarm_Jerk_Perc25 ),
                 data=input_data, method = "REML", algorithm = "gam")
summary(model_Perc25)

Family: gaussian 
Link function: identity 

Formula:
Y ~ c(xload) + c(xpace) + c(xgender) + c(xage) + c(xheight) + 
    c(xweight) + c(xhip) + c(xwaist) + +s(Wrist_Acc_Perc25) + 
    s(Wrist_Jerk_Perc25) + s(Trunk_Acc_Perc25) + s(Trunk_Jerk_Perc25) + 
    s(Upperarm_Acc_Perc25) + s(Upperarm_Jerk_Perc25)

Constant coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -2.223792   1.972056  -1.128 0.260295    
xload        0.100790   0.026937   3.742 0.000216 ***
xpace        0.872083   0.115794   7.531 4.95e-13 ***
xgenderF    -1.221753   0.448368  -2.725 0.006779 ** 
xage         0.083259   0.021461   3.880 0.000127 ***
xheight     -0.033624   0.012888  -2.609 0.009501 ** 
xweight      0.009688   0.007248   1.337 0.182276    
xhip         0.091905   0.009791   9.387  < 2e-16 ***
xwaist      -0.058944   0.012707  -4.639 5.08e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Smooth terms & functional coefficients:
                          edf Ref.df      F p-value    
Intercept(yindex)       7.301  9.000 59.588 < 2e-16 ***
s(Wrist_Acc_Perc25)     8.791 11.619  2.578 0.00301 ** 
s(Wrist_Jerk_Perc25)    1.000  1.000  0.016 0.90148    
s(Trunk_Acc_Perc25)     7.867 11.322  2.284 0.00992 ** 
s(Trunk_Jerk_Perc25)    2.090  2.549  3.305 0.02105 *  
s(Upperarm_Acc_Perc25)  1.875  2.313  5.962 0.00233 ** 
s(Upperarm_Jerk_Perc25) 5.306  6.166  3.372 0.00325 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.745   Deviance explained = 77.4%
-REML score = 454.87  Scale est. = 0.52451   n = 370 (37 x 10)

model_Perc75 <- pffr(Y ~
                 c(xload)  +
                 c(xpace)  +  
                 c(xgender)+ 
                 c(xage)+
                 c(xheight)+ 
                 c(xweight)+ 
                 c(xhip)+ 
                 c(xwaist)+
                 +s( Wrist_Acc_Perc75 )
                 +s( Wrist_Jerk_Perc75 )
                 +s( Trunk_Acc_Perc75 )
                 +s( Trunk_Jerk_Perc75 )
                 +s( Upperarm_Acc_Perc75 )
                 +s( Upperarm_Jerk_Perc75 ),
                 data=input_data, method = "REML", algorithm = "gam")
summary(model_Perc75)

Family: gaussian 
Link function: identity 

Formula:
Y ~ c(xload) + c(xpace) + c(xgender) + c(xage) + c(xheight) + 
    c(xweight) + c(xhip) + c(xwaist) + +s(Wrist_Acc_Perc75) + 
    s(Wrist_Jerk_Perc75) + s(Trunk_Acc_Perc75) + s(Trunk_Jerk_Perc75) + 
    s(Upperarm_Acc_Perc75) + s(Upperarm_Jerk_Perc75)

Constant coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.661929   2.333296   1.998   0.0466 *  
xload        0.032308   0.029752   1.086   0.2784    
xpace        0.769425   0.112412   6.845 4.29e-11 ***
xgenderF    -2.314791   0.482232  -4.800 2.50e-06 ***
xage         0.151736   0.024115   6.292 1.10e-09 ***
xheight     -0.088944   0.015502  -5.738 2.34e-08 ***
xweight      0.042807   0.007869   5.440 1.10e-07 ***
xhip         0.085847   0.009605   8.938  < 2e-16 ***
xwaist      -0.054999   0.013008  -4.228 3.13e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Smooth terms & functional coefficients:
                           edf Ref.df      F p-value    
Intercept(yindex)        7.739  9.000 83.091 < 2e-16 ***
s(Wrist_Acc_Perc75)     10.941 15.177  1.893 0.02272 *  
s(Wrist_Jerk_Perc75)     4.713  5.571  2.913 0.02738 *  
s(Trunk_Acc_Perc75)     13.926 18.671  4.935 < 2e-16 ***
s(Trunk_Jerk_Perc75)     9.520 13.225  2.577 0.00203 ** 
s(Upperarm_Acc_Perc75)   2.634  3.703  0.545 0.67970    
s(Upperarm_Jerk_Perc75)  9.867 13.653  1.160 0.30983    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.802   Deviance explained = 83.8%
-REML score = 433.47  Scale est. = 0.40775   n = 370 (37 x 10)

model_min <- pffr(Y ~
                 c(xload)  +  
                 c(xpace)  +  
                 c(xgender)+ 
                 c(xage)+
                 c(xheight)+ 
                 c(xweight)+ 
                 c(xhip)+ 
                 c(xwaist)+
                 +s( Wrist_Acc_Min )
                 +s( Wrist_Jerk_Min )
                 +s( Trunk_Acc_Min )
                 +s( Trunk_Jerk_Min )
                 +s( Upperarm_Acc_Min )
                 +s( Upperarm_Jerk_Min ),
                 data=input_data, method = "REML", algorithm = "gam")
summary(model_min)

Family: gaussian 
Link function: identity 

Formula:
Y ~ c(xload) + c(xpace) + c(xgender) + c(xage) + c(xheight) + 
    c(xweight) + c(xhip) + c(xwaist) + +s(Wrist_Acc_Min) + s(Wrist_Jerk_Min) + 
    s(Trunk_Acc_Min) + s(Trunk_Jerk_Min) + s(Upperarm_Acc_Min) + 
    s(Upperarm_Jerk_Min)

Constant coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  5.258278   1.986916   2.646  0.00851 ** 
xload        0.094120   0.014843   6.341 7.30e-10 ***
xpace        0.746249   0.107988   6.910 2.41e-11 ***
xgenderF    -2.367525   0.427075  -5.544 5.97e-08 ***
xage         0.110460   0.019364   5.704 2.55e-08 ***
xheight     -0.067324   0.011932  -5.642 3.55e-08 ***
xweight      0.016408   0.005982   2.743  0.00641 ** 
xhip         0.090582   0.009299   9.741  < 2e-16 ***
xwaist      -0.082668   0.011981  -6.900 2.57e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Smooth terms & functional coefficients:
                       edf Ref.df      F p-value    
Intercept(yindex)    7.161  9.000 55.313  <2e-16 ***
s(Wrist_Acc_Min)     2.990  3.590 12.299  <2e-16 ***
s(Wrist_Jerk_Min)    1.000  1.001  0.001  0.9813    
s(Trunk_Acc_Min)     2.563  3.124  2.957  0.0310 *  
s(Trunk_Jerk_Min)    1.000  1.000  5.627  0.0182 *  
s(Upperarm_Acc_Min)  6.539  9.597  1.887  0.0479 *  
s(Upperarm_Jerk_Min) 1.308  1.558  1.548  0.1384    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.726   Deviance explained = 74.9%
-REML score = 454.42  Scale est. = 0.5633    n = 370 (37 x 10)

model_max <- pffr(Y ~
                 c(xload)  +
                 c(xpace)  +  
                 c(xgender)+ 
                 c(xage)+
                 c(xheight)+ 
                 c(xweight)+ 
                 c(xhip)+ 
                 c(xwaist)+
                 +s( Wrist_Acc_Max )
                 +s( Wrist_Jerk_Max )
                 +s( Trunk_Acc_Max )
                 +s( Trunk_Jerk_Max )
                 +s( Upperarm_Acc_Max )
                 +s( Upperarm_Jerk_Max ),
                 data=input_data, method = "REML", algorithm = "gam")
summary(model_max)

Family: gaussian 
Link function: identity 

Formula:
Y ~ c(xload) + c(xpace) + c(xgender) + c(xage) + c(xheight) + 
    c(xweight) + c(xhip) + c(xwaist) + +s(Wrist_Acc_Max) + s(Wrist_Jerk_Max) + 
    s(Trunk_Acc_Max) + s(Trunk_Jerk_Max) + s(Upperarm_Acc_Max) + 
    s(Upperarm_Jerk_Max)

Constant coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.795067   2.037673   1.372 0.171094    
xload        0.094982   0.011953   7.947 3.09e-14 ***
xpace        0.763840   0.108138   7.064 9.73e-12 ***
xgenderF    -1.799837   0.408628  -4.405 1.44e-05 ***
xage         0.103169   0.019282   5.351 1.65e-07 ***
xheight     -0.059392   0.012020  -4.941 1.24e-06 ***
xweight      0.021589   0.006323   3.414 0.000719 ***
xhip         0.087727   0.009133   9.606  < 2e-16 ***
xwaist      -0.071740   0.012045  -5.956 6.64e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Smooth terms & functional coefficients:
                       edf Ref.df      F p-value    
Intercept(yindex)    7.256  9.000 58.974  <2e-16 ***
s(Wrist_Acc_Max)     2.056  2.533  4.145  0.0142 *  
s(Wrist_Jerk_Max)    1.001  1.001  0.580  0.4470    
s(Trunk_Acc_Max)     5.673  7.126  2.165  0.0353 *  
s(Trunk_Jerk_Max)    8.117 11.549  2.089  0.0181 *  
s(Upperarm_Acc_Max)  2.780  3.325  3.271  0.0213 *  
s(Upperarm_Jerk_Max) 5.370  7.497  1.632  0.1373    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.736   Deviance explained = 76.5%
-REML score = 457.01  Scale est. = 0.54305   n = 370 (37 x 10)

model_mad <- pffr(Y ~
                 c(xload)  +
                 c(xpace)  +  
                 c(xgender)+ 
                 c(xage)+
                 c(xheight)+ 
                 c(xweight)+ 
                 c(xhip)+ 
                 c(xwaist)+
                 +s( Wrist_Acc_Mad )
                 +s( Wrist_Jerk_Mad )
                 +s( Trunk_Acc_Mad )
                 +s( Trunk_Jerk_Mad )
                 +s( Upperarm_Acc_Mad )
                 +s( Upperarm_Jerk_Mad ),
                 data=input_data, method = "REML", algorithm = "gam")
summary(model_mad)

Family: gaussian 
Link function: identity 

Formula:
Y ~ c(xload) + c(xpace) + c(xgender) + c(xage) + c(xheight) + 
    c(xweight) + c(xhip) + c(xwaist) + +s(Wrist_Acc_Mad) + s(Wrist_Jerk_Mad) + 
    s(Trunk_Acc_Mad) + s(Trunk_Jerk_Mad) + s(Upperarm_Acc_Mad) + 
    s(Upperarm_Jerk_Mad)

Constant coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.281061   2.040259   1.118 0.264388    
xload       -0.025017   0.027394  -0.913 0.361819    
xpace        0.752193   0.113266   6.641 1.32e-10 ***
xgenderF    -1.289188   0.454835  -2.834 0.004881 ** 
xage         0.063662   0.020796   3.061 0.002389 ** 
xheight     -0.042516   0.012660  -3.358 0.000878 ***
xweight      0.016796   0.005972   2.813 0.005217 ** 
xhip         0.083488   0.009657   8.646 2.54e-16 ***
xwaist      -0.064180   0.012033  -5.334 1.82e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Smooth terms & functional coefficients:
                        edf Ref.df      F  p-value    
Intercept(yindex)     7.422  9.000 66.399  < 2e-16 ***
s(Wrist_Acc_Mad)     15.125 20.141  3.097 1.35e-05 ***
s(Wrist_Jerk_Mad)     1.000  1.000  2.460  0.11775    
s(Trunk_Acc_Mad)      1.353  1.608  0.555  0.61833    
s(Trunk_Jerk_Mad)     6.748  9.717  2.735  0.00345 ** 
s(Upperarm_Acc_Mad)   4.023  5.151  1.545  0.17959    
s(Upperarm_Jerk_Mad)  2.729  3.249  1.915  0.10837    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.758   Deviance explained = 78.9%
-REML score = 446.76  Scale est. = 0.49678   n = 370 (37 x 10)

model_rms <- pffr(Y ~
                 c(xload)  +
                 c(xpace)  +  
                 c(xgender)+ 
                 c(xage)+
                 c(xheight)+ 
                 c(xweight)+ 
                 c(xhip)+ 
                 c(xwaist)+
                 +s( Wrist_Acc_Rms )
                 +s( Wrist_Jerk_Rms )
                 +s( Trunk_Acc_Rms )
                 +s( Trunk_Jerk_Rms )
                 +s( Upperarm_Acc_Rms )
                 +s( Upperarm_Jerk_Rms ),
                 data=input_data, method = "REML", algorithm = "gam")
summary(model_rms)

Family: gaussian 
Link function: identity 

Formula:
Y ~ c(xload) + c(xpace) + c(xgender) + c(xage) + c(xheight) + 
    c(xweight) + c(xhip) + c(xwaist) + +s(Wrist_Acc_Rms) + s(Wrist_Jerk_Rms) + 
    s(Trunk_Acc_Rms) + s(Trunk_Jerk_Rms) + s(Upperarm_Acc_Rms) + 
    s(Upperarm_Jerk_Rms)

Constant coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  6.140008   2.065393   2.973  0.00320 ** 
xload        0.055843   0.017918   3.117  0.00201 ** 
xpace        0.966035   0.102345   9.439  < 2e-16 ***
xgenderF    -2.535610   0.398732  -6.359 7.75e-10 ***
xage         0.136807   0.019741   6.930 2.67e-11 ***
xheight     -0.091324   0.012498  -7.307 2.59e-12 ***
xweight      0.039107   0.006912   5.658 3.64e-08 ***
xhip         0.101164   0.009241  10.947  < 2e-16 ***
xwaist      -0.085127   0.012129  -7.018 1.56e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Smooth terms & functional coefficients:
                        edf Ref.df      F  p-value    
Intercept(yindex)     7.800  9.000 91.504  < 2e-16 ***
s(Wrist_Acc_Rms)     15.096 19.214  3.888  < 2e-16 ***
s(Wrist_Jerk_Rms)     6.546  7.837  4.738 2.82e-05 ***
s(Trunk_Acc_Rms)     10.641 14.463  5.601  < 2e-16 ***
s(Trunk_Jerk_Rms)     8.998 12.686  2.465 0.003696 ** 
s(Upperarm_Acc_Rms)  14.755 19.153  2.778 0.000124 ***
s(Upperarm_Jerk_Rms)  3.507  4.010  2.873 0.022298 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.838   Deviance explained = 87.1%
-REML score = 412.26  Scale est. = 0.33204   n = 370 (37 x 10)

model_std <- pffr(Y ~
                 c(xload)  +
                 c(xpace)  +  
                 c(xgender)+ 
                 c(xage)+
                 c(xheight)+ 
                 c(xweight)+ 
                 c(xhip)+ 
                 c(xwaist)+
                 +s( Wrist_Acc_Std )
                 +s( Wrist_Jerk_Std )
                 +s( Trunk_Acc_Std )
                 +s( Trunk_Jerk_Std )
                 +s( Upperarm_Acc_Std )
                 +s( Upperarm_Jerk_Std ),
                 data=input_data, method = "REML", algorithm = "gam")
summary(model_std)

Family: gaussian 
Link function: identity 

Formula:
Y ~ c(xload) + c(xpace) + c(xgender) + c(xage) + c(xheight) + 
    c(xweight) + c(xhip) + c(xwaist) + +s(Wrist_Acc_Std) + s(Wrist_Jerk_Std) + 
    s(Trunk_Acc_Std) + s(Trunk_Jerk_Std) + s(Upperarm_Acc_Std) + 
    s(Upperarm_Jerk_Std)

Constant coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  5.511079   2.123564   2.595 0.009884 ** 
xload        0.020105   0.018080   1.112 0.266958    
xpace        0.804359   0.110922   7.252 3.06e-12 ***
xgenderF    -2.135145   0.444383  -4.805 2.38e-06 ***
xage         0.103567   0.019017   5.446 1.02e-07 ***
xheight     -0.068319   0.012747  -5.360 1.59e-07 ***
xweight      0.023340   0.006098   3.827 0.000156 ***
xhip         0.095626   0.008885  10.763  < 2e-16 ***
xwaist      -0.086294   0.011719  -7.364 1.50e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Smooth terms & functional coefficients:
                       edf Ref.df      F p-value    
Intercept(yindex)    7.544  9.000 69.737 < 2e-16 ***
s(Wrist_Acc_Std)     3.569  4.142 12.410 < 2e-16 ***
s(Wrist_Jerk_Std)    8.086 11.612  1.420 0.15286    
s(Trunk_Acc_Std)     5.359  6.776  1.609 0.13820    
s(Trunk_Jerk_Std)    7.831 11.158  2.543 0.00422 ** 
s(Upperarm_Acc_Std)  3.014  3.892  0.883 0.39434    
s(Upperarm_Jerk_Std) 2.238  2.700  1.610 0.19463    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.771   Deviance explained =   80%
-REML score = 436.15  Scale est. = 0.46998   n = 370 (37 x 10)

3.5 Modeling with raw sensors as functional inputs

In this section we extracted functional features from the raw sensor signals. In order to do that, we first applied a 4th order Butterworth filter in the signals and then the filtered signals were used to extract the features.

3.5.1 Data Extraction

y_T_val = t(y_T_val)

Conditions = as.vector(t(design[4:7, seq(5,69,4)]))
ActiGraph_IMU_DF$Condition = rep(Conditions, 3)

num_sub = length(unique(ActiGraph_IMU_DF$Subject))
num_session = length(unique(ActiGraph_IMU_DF$Session))

Session_List = vector(mode = "list", length = num_session)

wrist_total_filtered = vector(mode = "list", length = num_sub)
trunk_total_filtered = vector(mode = "list", length = num_sub)
upperarm_total_filtered = vector(mode = "list", length = num_sub)

wrist_total = vector(mode = "list", length = num_sub)
trunk_total = vector(mode = "list", length = num_sub)
upperarm_total = vector(mode = "list", length = num_sub)

for (i in 1:num_sub) {
  wrist_total_filtered[[i]] = Session_List
  trunk_total_filtered[[i]] = Session_List
  upperarm_total_filtered[[i]] = Session_List
  wrist_total[[i]] = Session_List
  trunk_total[[i]] = Session_List
  upperarm_total[[i]] = Session_List
  
}

# Butterworth Filter Function
butter_filter_data = function(data){
  n = 4 # The order of the Butterworth filter, which determines the complexity of the filter's transfer function.
  fc = 0.2 # The cutoff frequency which is the frequency at which the filter's response starts to attenuate.
  Fs = 100 # The sampling rate is 100 Hz for all IMU signals.
  Wn = (2/Fs)*fc
  butter_filter = signal::butter(n=n, W=Wn, type="low")
  filt_data = signal::filtfilt(data, filt = butter_filter)
}

# Data frame filter function

df_filter = function(df, threshold=15, metric = 'median', sub, id){
  metric = sym(metric)
  df |>
    dplyr::mutate_at(2:10,.funs=butter_filter_data) |>
    dplyr::mutate(time= as.POSIXct(df[,1], origin="1970-01-01")) |>
    dplyr::mutate(time_diff_secs= as.numeric(difftime(time, time[1], units = "secs"))) |>
    dplyr::mutate(time_diff_mins= round(as.numeric(difftime(time, time[1], units = "mins")), digits = 0)) |>
    dplyr::mutate(max_time_diff_secs = max(time_diff_secs, na.rm = TRUE)) |>
    dplyr::mutate(interval = round(time_diff_secs / 300, digits = 0) * 300) |>
    dplyr::filter(abs(time_diff_secs - interval) <= threshold | (max_time_diff_secs-time_diff_secs)<=(2*threshold), time_diff_mins<=45)|>
    dplyr::relocate(c(time,time_diff_secs,time_diff_mins), .after = Timestamp)|>
    dplyr::mutate(interval = round(time_diff_mins / 5, digits = 0)) |>
    dplyr::group_by(interval) |>
    dplyr::summarise(across(1:13, !!metric, na.rm = TRUE)) |>
    dplyr::ungroup() |>
    dplyr::mutate(time = time_diff_mins, subject = sub, condition = substr(ActiGraph_IMU_DF$Condition [id],11,16), sub_con = paste(subject,condition, sep="-")) |>
    dplyr::relocate(c(subject,condition,sub_con), .after = time)|>
    dplyr::select(-interval,-time_diff_mins,-time_diff_secs)|>
    mutate(time = if_else(row_number() == n() & last(time) != 45, 45, time))
}

for (sub in 1:num_sub) {
  for (session in 1:num_session) {
  if (!is.na(CP$CP1[(session-1)*num_sub+sub]))
  {
  tryCatch({
    
  trunk_id = which(ActiGraph_IMU_DF$Subject == paste0("Sub", str_pad(sub, 2, pad = "0")) & 
  ActiGraph_IMU_DF$Session == paste0("Session", session) & ActiGraph_IMU_DF$Part == "trunk")

  wrist_id = which(ActiGraph_IMU_DF$Subject == paste0("Sub", str_pad(sub, 2, pad = "0")) & 
                   ActiGraph_IMU_DF$Session == paste0("Session", session) & ActiGraph_IMU_DF$Part == "wrist")
  
  upperarm_id = which(ActiGraph_IMU_DF$Subject == paste0("Sub", str_pad(sub, 2, pad = "0")) & 
                   ActiGraph_IMU_DF$Session == paste0("Session", session) & ActiGraph_IMU_DF$Part == "upper arm")

  CP1 = CP[CP$Subject == paste0("Sub", str_pad(sub, 2, pad = "0")) & CP$Session == paste0("Session", session), 3]
  CP2 = CP[CP$Subject == paste0("Sub", str_pad(sub, 2, pad = "0")) & CP$Session == paste0("Session", session), 4]
  

  trunk_data <- ActiGraph_IMU_Data[[trunk_id]]
  trunk_data <- trunk_data[which(trunk_data$Timestamp >= CP1 & trunk_data$Timestamp <= CP2),]
  
  upperarm_data <- ActiGraph_IMU_Data[[upperarm_id]]
  upperarm_data = upperarm_data[which(upperarm_data$Timestamp >= CP1 & upperarm_data$Timestamp <= CP2),]
  
  wrist_data <- ActiGraph_IMU_Data[[wrist_id]]
  wrist_data = wrist_data[which(wrist_data$Timestamp >= CP1 & wrist_data$Timestamp <= CP2),]
  

  if ((sub == 15 & session == 3) | sub == 16) {
    trunk_data = trunk_data |> dplyr::mutate(trunk_data[,c(2,3,5,6,8,9)],-1*trunk_data[,c(2,3,5,6,8,9)])
  }
  
  if(is.null(trunk_data)) {
  trunk_data_filtered <- NULL
  trunk_total_filtered[[sub]][[session]]<- NULL
  trunk_total[[sub]][[session]] <- NULL
  } else {
    trunk_total[[sub]][[session]] = trunk_data |> dplyr::mutate(sub_con = paste(sub,substr(ActiGraph_IMU_DF$Condition [trunk_id],11,16),sep = "-")) |>
      dplyr::relocate(sub_con,.after = Timestamp)
                                                                 
    trunk_data_filtered <-df_filter(trunk_data, threshold=15, metric = 'median', sub, trunk_id)
    trunk_total_filtered[[sub]][[session]] = trunk_data_filtered
    names(trunk_total_filtered[[sub]])[session] <- substr(ActiGraph_IMU_DF$Condition [trunk_id],11,16)
  }
  
  if(is.null(wrist_data)) {
  wrist_data_filtered <- NULL
  wrist_total_filtered[[sub]][[session]] <- NULL
  wrist_total[[sub]][[session]] <- NULL
  } else {
    wrist_total[[sub]][[session]] = wrist_data  |>  dplyr::mutate(sub_con = paste(sub,substr(ActiGraph_IMU_DF$Condition [wrist_id],11,16),sep = "-")) |>
      dplyr::relocate(sub_con,.after = Timestamp)
    
    wrist_data_filtered <-df_filter(wrist_data, threshold=15, metric = 'median', sub, wrist_id)
    wrist_total_filtered[[sub]][[session]] = wrist_data_filtered
    names(wrist_total_filtered[[sub]])[session] <- substr(ActiGraph_IMU_DF$Condition[wrist_id],11,16)
  }
  
  if(is.null(upperarm_data)) {
  upperarm_data_filtered <- NULL
  upperarm_total_filtered[[sub]][[session]] <- NULL
  upperarm_total[[sub]][[session]] <- NULL
  } else {
    upperarm_total[[sub]][[session]] = upperarm_data |> dplyr::mutate(sub_con = paste(sub,substr(ActiGraph_IMU_DF$Condition [upperarm_id],11,16),sep = "-")) |>
      dplyr::relocate(sub_con,.after = Timestamp)
    
    upperarm_data_filtered = df_filter(upperarm_data, threshold=15, metric = 'median', sub, upperarm_id)
    upperarm_total_filtered[[sub]][[session]] = upperarm_data_filtered
    names(upperarm_total_filtered[[sub]])[session] <- substr(ActiGraph_IMU_DF$Condition[upperarm_id],11,16)
  }
  
  }, error = function(e) {print(e)})
  }
  }
}

trunk_total_df = lapply(trunk_total, dplyr::bind_rows)
trunk_total_df = dplyr::bind_rows(trunk_total_df)
trunk_total_df = trunk_total_df |> 
  dplyr::filter(sub_con %in% colnames(y_T_val)) |> 
  dplyr::rename_with(~paste0('trunk_',.), .cols = 3:11) |>
  dplyr::rename_with(~gsub("\\.", "_", .), .cols = 3:11)

trunk_total_filtered_df = lapply(trunk_total_filtered, dplyr::bind_rows)
trunk_total_filtered_df = dplyr::bind_rows(trunk_total_filtered_df)
trunk_total_filtered_df = trunk_total_filtered_df |> 
  dplyr::filter(sub_con %in% colnames(y_T_val)) |> 
  dplyr::rename_with(~paste0('trunk_',.), .cols = 6:14) |>
  dplyr::rename_with(~gsub("\\.", "_", .), .cols = 6:14)

wrist_total_df = lapply(wrist_total, dplyr::bind_rows)
wrist_total_df = dplyr::bind_rows(wrist_total_df)
wrist_total_df = wrist_total_df |> 
  dplyr::filter(sub_con %in% colnames(y_T_val)) |> 
  dplyr::rename_with(~paste0('wrist_',.), .cols = 3:11) |>
  dplyr::rename_with(~gsub("\\.", "_", .), .cols = 3:11)

wrist_total_filtered_df = lapply(wrist_total_filtered, dplyr::bind_rows)
wrist_total_filtered_df = dplyr::bind_rows(wrist_total_filtered_df)
wrist_total_filtered_df = wrist_total_filtered_df |> 
  dplyr::filter(sub_con %in% colnames(y_T_val)) |> 
  dplyr::rename_with(~paste0('wrist_',.), .cols = 6:14) |>
  dplyr::rename_with(~gsub("\\.", "_", .), .cols = 6:14)

upperarm_total_df = lapply(upperarm_total, dplyr::bind_rows)
upperarm_total_df = dplyr::bind_rows(upperarm_total_df)
upperarm_total_df = upperarm_total_df |> 
  dplyr::filter(sub_con %in% colnames(y_T_val)) |> 
  dplyr::rename_with(~paste0('upperarm_',.), .cols = 3:11) |>
  dplyr::rename_with(~gsub("\\.", "_", .), .cols = 3:11)

upperarm_total_filtered_df = lapply(upperarm_total_filtered, dplyr::bind_rows)
upperarm_total_filtered_df = dplyr::bind_rows(upperarm_total_filtered_df)
upperarm_total_filtered_df = upperarm_total_filtered_df |> 
  dplyr::filter(sub_con %in% colnames(y_T_val)) |>
  dplyr::rename_with(~paste0('upperarm_',.), .cols = 6:14) |>
  dplyr::rename_with(~gsub("\\.", "_", .), .cols = 6:14)

rm(trunk_data_filtered,trunk_total,trunk_total_filtered, wrist_total,wrist_total_filtered, upperarm_total, upperarm_total_filtered )

# sensors_total_df = trunk_total_df |>
#   dplyr::left_join(upperarm_total_df, by = c('sub_con','Timestamp')) |>
#   dplyr::left_join(wrist_total_filtered_df, by = c('sub_con','Timestamp'))
# Got this error: Error: vector memory exhausted (limit reached?)

sensors_total_filtered_df = trunk_total_filtered_df |>
  dplyr::left_join(upperarm_total_filtered_df, by = c('sub_con','time','subject','condition','Timestamp')) |>
  dplyr::left_join(wrist_total_filtered_df, by = c('sub_con','time','subject','condition','Timestamp'))

save(list = c("trunk_total_df", "trunk_total_filtered_df", "wrist_total_df", "wrist_total_filtered_df", "upperarm_total_df", "upperarm_total_filtered_df", "sensors_total_filtered_df"), file = "sensors.RData")

subjects_sensors = substr(colnames(y_T_val),1,2)
subjects_sensors[which(substr(subjects_sensors,2,2)=="-")]=paste0(0,substr(subjects_sensors[which(substr(subjects_sensors,2,2)=="-")],1,1))
subjects_sensors = paste0("Sub", subjects_sensors)

input_data_sensors = vector(mode= "list")
input_data_sensors[[1]] <- input_data_sensors$Y <- t(y_T_val)
input_data_sensors[[2]] <- input_data_sensors$xload <- c(rep(2.5,31),rep(1.5,13))
input_data_sensors[[3]] <- input_data_sensors$xpace <- c(rep(5,11),rep(10,10),rep(15,23))
input_data_sensors[[4]] <- input_data_sensors$xgender <- factor(Anthrpmtrc_data$Gender[match(subjects_sensors, Anthrpmtrc_data$`Sub#`)], levels=c("F","M"))
input_data_sensors[[5]] <- input_data_sensors$xage <- Anthrpmtrc_data$Age[match(subjects_sensors,Anthrpmtrc_data$`Sub#`)]
input_data_sensors[[6]] <- input_data_sensors$xheight <- Anthrpmtrc_data$`Height (cm)`[match(subjects_sensors,Anthrpmtrc_data$`Sub#`)]
input_data_sensors[[7]] <- input_data_sensors$xweight <- Anthrpmtrc_data$`Weight (kg)`[match(subjects_sensors,Anthrpmtrc_data$`Sub#`)]
input_data_sensors[[8]] <- input_data_sensors$xhip <- Anthrpmtrc_data$`Hip circumference (cm)`[match(subjects_sensors,Anthrpmtrc_data$`Sub#`)]
input_data_sensors[[9]] <- input_data_sensors$xwaist <- Anthrpmtrc_data$`Waist circumference (cm)`[match(subjects_sensors,Anthrpmtrc_data$`Sub#`)]

for (i in 1:27){
  name = colnames(sensors_total_filtered_df)[(5+i)]
  df_temp <- sensors_total_filtered_df[,c(2,5,(5+i))] 
  df_temp <- df_temp |> pivot_wider(names_from = time, values_from = name)
  df_temp<- df_temp[order(match(df_temp$sub_con, colnames(y_T_val))),]
  df_temp = df_temp[,-1]
  rownames(df_temp) <- colnames(y_T_val)
  df_temp = as.matrix(df_temp)
  
  input_data_sensors[[(9+i)]] <-  input_data_sensors[[name]] <- df_temp
}

3.5.2 Modeling

# Below model considers smoothed sensors data in addition to the individual characteristics and task specifications.

model_trunk_sensors <- pffr(Y ~
                 c(xload)  +
                 c(xpace)  +  
                 c(xgender)+ 
                 c(xage)+
                 c(xheight)+ 
                 c(xweight)+ 
                 c(xhip)+ 
                 c(xwaist)+
                 +s( trunk_Accelerometer_X )
                 +s( trunk_Accelerometer_Y )
                 +s( trunk_Accelerometer_Z )
                 +s( trunk_Gyroscope_X )
                 +s( trunk_Gyroscope_Y )
                 +s( trunk_Gyroscope_Z ),
                 # +s( trunk_Magnetometer_X )
                 # +s( trunk_Magnetometer_Y )
                 # +s( trunk_Magnetometer_Z ),
                 data=input_data_sensors, method = "REML", algorithm = "gam")

summary(model_trunk_sensors)

Family: gaussian 
Link function: identity 

Formula:
Y ~ c(xload) + c(xpace) + c(xgender) + c(xage) + c(xheight) + 
    c(xweight) + c(xhip) + c(xwaist) + +s(trunk_Accelerometer_X) + 
    s(trunk_Accelerometer_Y) + s(trunk_Accelerometer_Z) + s(trunk_Gyroscope_X) + 
    s(trunk_Gyroscope_Y) + s(trunk_Gyroscope_Z)

Constant coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -10.278754   1.106610  -9.289  < 2e-16 ***
xload         0.573909   0.106088   5.410 1.10e-07 ***
xpace         0.092282   0.012196   7.567 2.78e-13 ***
xgenderM     -1.748958   0.241387  -7.245 2.32e-12 ***
xage          0.042326   0.014284   2.963  0.00323 ** 
xheight       0.005630   0.008431   0.668  0.50467    
xweight       0.011452   0.006311   1.815  0.07035 .  
xhip          0.040027   0.006870   5.826 1.19e-08 ***
xwaist        0.045755   0.006342   7.215 2.83e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Smooth terms & functional coefficients:
                            edf Ref.df      F  p-value    
Intercept(yindex)         7.770  9.000 82.570  < 2e-16 ***
s(trunk_Accelerometer_X)  5.909  6.952  6.111 1.45e-06 ***
s(trunk_Accelerometer_Y)  2.640  3.748  0.856  0.47783    
s(trunk_Accelerometer_Z) 12.153 16.625  7.541  < 2e-16 ***
s(trunk_Gyroscope_X)      4.281  6.305  0.709  0.64439    
s(trunk_Gyroscope_Y)      1.000  1.001  8.109  0.00463 ** 
s(trunk_Gyroscope_Z)      6.475  9.369  1.358  0.20042    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.767   Deviance explained = 79.3%
-REML score = 547.77  Scale est. = 0.54849   n = 440 (44 x 10)

model_upperarm_sensors <- pffr(Y ~
                 c(xload)  +
                 c(xpace)  +  
                 c(xgender)+ 
                 c(xage)+
                 c(xheight)+ 
                 c(xweight)+ 
                 c(xhip)+ 
                 c(xwaist)+
                 +s( upperarm_Accelerometer_X )
                 +s( upperarm_Accelerometer_Y )
                 +s( upperarm_Accelerometer_Z )
                 +s( upperarm_Gyroscope_X )
                 +s( upperarm_Gyroscope_Y )
                 +s( upperarm_Gyroscope_Z ),
                 # +s( upperarm_Magnetometer_X )
                 # +s( upperarm_Magnetometer_Y )
                 # +s( upperarm_Magnetometer_Z ),
                 data=input_data_sensors, method = "REML", algorithm = "gam")

summary(model_upperarm_sensors)

Family: gaussian 
Link function: identity 

Formula:
Y ~ c(xload) + c(xpace) + c(xgender) + c(xage) + c(xheight) + 
    c(xweight) + c(xhip) + c(xwaist) + +s(upperarm_Accelerometer_X) + 
    s(upperarm_Accelerometer_Y) + s(upperarm_Accelerometer_Z) + 
    s(upperarm_Gyroscope_X) + s(upperarm_Gyroscope_Y) + s(upperarm_Gyroscope_Z)

Constant coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -9.598721   1.065387  -9.010  < 2e-16 ***
xload        0.764778   0.112647   6.789 4.41e-11 ***
xpace        0.110335   0.016078   6.863 2.80e-11 ***
xgenderM    -1.489485   0.280056  -5.319 1.80e-07 ***
xage         0.031732   0.015632   2.030 0.043072 *  
xheight      0.024496   0.007383   3.318 0.000995 ***
xweight      0.013478   0.006215   2.168 0.030748 *  
xhip         0.016572   0.007471   2.218 0.027135 *  
xwaist       0.020179   0.006560   3.076 0.002251 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Smooth terms & functional coefficients:
                               edf Ref.df      F  p-value    
Intercept(yindex)            7.529  9.000 67.070  < 2e-16 ***
s(upperarm_Accelerometer_X) 11.457 15.385  2.668 0.000641 ***
s(upperarm_Accelerometer_Y)  4.643  5.630  7.429 8.34e-07 ***
s(upperarm_Accelerometer_Z) 11.563 15.090  5.322  < 2e-16 ***
s(upperarm_Gyroscope_X)      5.022  6.771  1.505 0.164747    
s(upperarm_Gyroscope_Y)      5.119  7.186  1.480 0.188987    
s(upperarm_Gyroscope_Z)      9.367 13.093  1.922 0.026849 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.761   Deviance explained = 79.5%
-REML score = 569.63  Scale est. = 0.56274   n = 440 (44 x 10)

model_wrist_sensors <- pffr(Y ~
                 c(xload)  +
                 c(xpace)  +  
                 c(xgender)+ 
                 c(xage)+
                 c(xheight)+ 
                 c(xweight)+ 
                 c(xhip)+ 
                 c(xwaist)+
                 +s( wrist_Accelerometer_X )
                 +s( wrist_Accelerometer_Y )
                 +s( wrist_Accelerometer_Z )
                 +s( wrist_Gyroscope_X )
                 +s( wrist_Gyroscope_Y )
                 +s( wrist_Gyroscope_Z ),
                 # +s( wrist_Magnetometer_X )
                 # +s( wrist_Magnetometer_Y )
                 # +s( wrist_Magnetometer_Z ),
                 data=input_data_sensors, method = "REML", algorithm = "gam")

summary(model_wrist_sensors)

Family: gaussian 
Link function: identity 

Formula:
Y ~ c(xload) + c(xpace) + c(xgender) + c(xage) + c(xheight) + 
    c(xweight) + c(xhip) + c(xwaist) + +s(wrist_Accelerometer_X) + 
    s(wrist_Accelerometer_Y) + s(wrist_Accelerometer_Z) + s(wrist_Gyroscope_X) + 
    s(wrist_Gyroscope_Y) + s(wrist_Gyroscope_Z)

Constant coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -9.326992   1.042150  -8.950  < 2e-16 ***
xload        0.601678   0.109204   5.510 6.73e-08 ***
xpace        0.090530   0.016490   5.490 7.46e-08 ***
xgenderM    -1.161328   0.266224  -4.362 1.67e-05 ***
xage        -0.028449   0.014257  -1.995   0.0467 *  
xheight      0.029574   0.007380   4.007 7.42e-05 ***
xweight     -0.009198   0.006667  -1.380   0.1685    
xhip         0.052552   0.008125   6.468 3.13e-10 ***
xwaist       0.007742   0.007240   1.069   0.2856    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Smooth terms & functional coefficients:
                            edf Ref.df      F  p-value    
Intercept(yindex)         7.586  9.000 71.676  < 2e-16 ***
s(wrist_Accelerometer_X)  7.060  9.574  5.569 2.02e-07 ***
s(wrist_Accelerometer_Y)  6.930  9.736  3.410 0.000347 ***
s(wrist_Accelerometer_Z)  6.687  9.551  1.561 0.133440    
s(wrist_Gyroscope_X)     13.023 17.250  2.180 0.004599 ** 
s(wrist_Gyroscope_Y)      4.185  5.860  1.491 0.175779    
s(wrist_Gyroscope_Z)      3.817  5.579  1.233 0.297355    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =   0.77   Deviance explained = 80.1%
-REML score = 542.78  Scale est. = 0.54625   n = 430 (44 x 10)

model_acc_sensors <- pffr(Y ~
                 c(xload)  +
                 c(xpace)  +  
                 c(xgender)+ 
                 c(xage)+
                 c(xheight)+ 
                 c(xweight)+ 
                 c(xhip)+ 
                 c(xwaist)+
                 +s( trunk_Accelerometer_X )
                 +s( trunk_Accelerometer_Y )
                 +s( trunk_Accelerometer_Z )
                 +s( upperarm_Accelerometer_X )
                 +s( upperarm_Accelerometer_Y )
                 +s( upperarm_Accelerometer_Z )
                 +s( wrist_Accelerometer_X )
                 +s( wrist_Accelerometer_Y )
                 +s( wrist_Accelerometer_Z ),
                 data=input_data_sensors, method = "REML", algorithm = "gam")

summary(model_acc_sensors)

Family: gaussian 
Link function: identity 

Formula:
Y ~ c(xload) + c(xpace) + c(xgender) + c(xage) + c(xheight) + 
    c(xweight) + c(xhip) + c(xwaist) + +s(trunk_Accelerometer_X) + 
    s(trunk_Accelerometer_Y) + s(trunk_Accelerometer_Z) + s(upperarm_Accelerometer_X) + 
    s(upperarm_Accelerometer_Y) + s(upperarm_Accelerometer_Z) + 
    s(wrist_Accelerometer_X) + s(wrist_Accelerometer_Y) + s(wrist_Accelerometer_Z)

Constant coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -9.574819   0.863508 -11.088  < 2e-16 ***
xload        0.958402   0.081092  11.819  < 2e-16 ***
xpace        0.107272   0.013322   8.052 1.69e-14 ***
xgenderM    -1.660544   0.209122  -7.941 3.58e-14 ***
xage         0.078385   0.013318   5.886 1.01e-08 ***
xheight     -0.002126   0.007180  -0.296 0.767316    
xweight      0.020721   0.005807   3.568 0.000416 ***
xhip         0.031686   0.005742   5.519 7.14e-08 ***
xwaist       0.031368   0.005778   5.429 1.13e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Smooth terms & functional coefficients:
                               edf Ref.df       F  p-value    
Intercept(yindex)            8.340  9.000 152.657  < 2e-16 ***
s(trunk_Accelerometer_X)     9.079 11.330   3.986 1.72e-05 ***
s(trunk_Accelerometer_Y)    13.575 16.828   2.835 0.000195 ***
s(trunk_Accelerometer_Z)    12.795 16.993   5.801  < 2e-16 ***
s(upperarm_Accelerometer_X)  9.421 12.385   1.729 0.059532 .  
s(upperarm_Accelerometer_Y) 12.548 15.276   7.265  < 2e-16 ***
s(upperarm_Accelerometer_Z) 16.258 20.673   5.544  < 2e-16 ***
s(wrist_Accelerometer_X)    13.916 18.126   6.757  < 2e-16 ***
s(wrist_Accelerometer_Y)     1.001  1.001  12.759 0.000408 ***
s(wrist_Accelerometer_Z)     9.022 11.432   4.408 2.75e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.912   Deviance explained = 93.6%
-REML score = 413.12  Scale est. = 0.20828   n = 430 (44 x 10)

#The model with only the functional inputs.
model_acc2_sensors <- pffr(Y ~
                 +s( trunk_Accelerometer_X )
                 +s( trunk_Accelerometer_Y )
                 +s( trunk_Accelerometer_Z )
                 +s( upperarm_Accelerometer_X )
                 +s( upperarm_Accelerometer_Y )
                 +s( upperarm_Accelerometer_Z )
                 +s( wrist_Accelerometer_X )
                 +s( wrist_Accelerometer_Y )
                 +s( wrist_Accelerometer_Z ),
                 data=input_data_sensors, method = "REML", algorithm = "gam")

summary(model_acc2_sensors)

Family: gaussian 
Link function: identity 

Formula:
Y ~ +s(trunk_Accelerometer_X) + s(trunk_Accelerometer_Y) + s(trunk_Accelerometer_Z) + 
    s(upperarm_Accelerometer_X) + s(upperarm_Accelerometer_Y) + 
    s(upperarm_Accelerometer_Z) + s(wrist_Accelerometer_X) + 
    s(wrist_Accelerometer_Y) + s(wrist_Accelerometer_Z)

Constant coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.91709    0.04309    44.5   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Smooth terms & functional coefficients:
                               edf Ref.df      F  p-value    
Intercept(yindex)            7.227  9.000 48.896  < 2e-16 ***
s(trunk_Accelerometer_X)     5.293  6.897  2.284 0.029515 *  
s(trunk_Accelerometer_Y)    10.080 12.643  2.893 0.000432 ***
s(trunk_Accelerometer_Z)     8.376 11.926  3.142 0.000272 ***
s(upperarm_Accelerometer_X)  2.126  2.577  0.918 0.364474    
s(upperarm_Accelerometer_Y)  5.234  6.027  3.797 0.000933 ***
s(upperarm_Accelerometer_Z) 13.528 17.171  4.225  < 2e-16 ***
s(wrist_Accelerometer_X)     8.892 12.382  5.068  < 2e-16 ***
s(wrist_Accelerometer_Y)     2.523  3.013  8.797 1.22e-05 ***
s(wrist_Accelerometer_Z)     4.976  7.029  1.842 0.078686 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.683   Deviance explained = 73.3%
-REML score = 605.28  Scale est. = 0.75247   n = 430 (44 x 10)

model_gyr_sensors <- pffr(Y ~
                 c(xload)  +
                 c(xpace)  +  
                 c(xgender)+ 
                 c(xage)+
                 c(xheight)+ 
                 c(xweight)+ 
                 c(xhip)+ 
                 c(xwaist)+
                 +s( trunk_Gyroscope_X )
                 +s( trunk_Gyroscope_Y )
                 +s( trunk_Gyroscope_Z )
                 +s( upperarm_Gyroscope_X )
                 +s( upperarm_Gyroscope_Y )
                 +s( upperarm_Gyroscope_Z )
                 +s( wrist_Gyroscope_X )
                 +s( wrist_Gyroscope_Y )
                 +s( wrist_Gyroscope_Z ),
                 data=input_data_sensors, method = "REML", algorithm = "gam")

summary(model_gyr_sensors)

Family: gaussian 
Link function: identity 

Formula:
Y ~ c(xload) + c(xpace) + c(xgender) + c(xage) + c(xheight) + 
    c(xweight) + c(xhip) + c(xwaist) + +s(trunk_Gyroscope_X) + 
    s(trunk_Gyroscope_Y) + s(trunk_Gyroscope_Z) + s(upperarm_Gyroscope_X) + 
    s(upperarm_Gyroscope_Y) + s(upperarm_Gyroscope_Z) + s(wrist_Gyroscope_X) + 
    s(wrist_Gyroscope_Y) + s(wrist_Gyroscope_Z)

Constant coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -8.255193   1.101463  -7.495 4.79e-13 ***
xload        0.527681   0.118909   4.438 1.20e-05 ***
xpace        0.100326   0.013179   7.613 2.18e-13 ***
xgenderM    -1.225559   0.267254  -4.586 6.17e-06 ***
xage         0.006938   0.014747   0.471   0.6383    
xheight      0.015561   0.007953   1.957   0.0511 .  
xweight      0.012844   0.006859   1.873   0.0619 .  
xhip         0.041325   0.008234   5.019 8.02e-07 ***
xwaist       0.007540   0.006981   1.080   0.2808    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Smooth terms & functional coefficients:
                          edf Ref.df      F  p-value    
Intercept(yindex)       7.250  9.000 62.467  < 2e-16 ***
s(trunk_Gyroscope_X)    7.570 10.704  1.384   0.1847    
s(trunk_Gyroscope_Y)    5.689  8.269  1.862   0.0561 .  
s(trunk_Gyroscope_Z)    1.000  1.000  4.517   0.0342 *  
s(upperarm_Gyroscope_X) 5.113  7.096  0.909   0.4926    
s(upperarm_Gyroscope_Y) 5.366  7.165  4.970 1.98e-05 ***
s(upperarm_Gyroscope_Z) 2.716  3.943  0.619   0.6388    
s(wrist_Gyroscope_X)    6.137  8.433  1.989   0.0437 *  
s(wrist_Gyroscope_Y)    2.608  3.327  2.838   0.0324 *  
s(wrist_Gyroscope_Z)    1.001  1.002  5.234   0.0226 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.729   Deviance explained = 76.2%
-REML score = 564.86  Scale est. = 0.64336   n = 430 (44 x 10)

#The model with only the functional inputs.
model_gyr2_sensors <- pffr(Y ~
                 +s( trunk_Gyroscope_X )
                 +s( trunk_Gyroscope_Y )
                 +s( trunk_Gyroscope_Z )
                 +s( upperarm_Gyroscope_X )
                 +s( upperarm_Gyroscope_Y )
                 +s( upperarm_Gyroscope_Z )
                 +s( wrist_Gyroscope_X )
                 +s( wrist_Gyroscope_Y )
                 +s( wrist_Gyroscope_Z ),
                 data=input_data_sensors, method = "REML", algorithm = "gam")

summary(model_gyr2_sensors)

Family: gaussian 
Link function: identity 

Formula:
Y ~ +s(trunk_Gyroscope_X) + s(trunk_Gyroscope_Y) + s(trunk_Gyroscope_Z) + 
    s(upperarm_Gyroscope_X) + s(upperarm_Gyroscope_Y) + s(upperarm_Gyroscope_Z) + 
    s(wrist_Gyroscope_X) + s(wrist_Gyroscope_Y) + s(wrist_Gyroscope_Z)

Constant coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.94204    0.05166   37.59   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Smooth terms & functional coefficients:
                          edf Ref.df      F p-value    
Intercept(yindex)       6.757  9.000 39.174  <2e-16 ***
s(trunk_Gyroscope_X)    6.684  9.579  1.204  0.2929    
s(trunk_Gyroscope_Y)    1.911  2.332  1.254  0.2565    
s(trunk_Gyroscope_Z)    3.183  4.416  0.789  0.5457    
s(upperarm_Gyroscope_X) 1.551  1.924  0.954  0.3178    
s(upperarm_Gyroscope_Y) 4.891  6.267  1.789  0.0950 .  
s(upperarm_Gyroscope_Z) 8.496 11.814  1.368  0.1939    
s(wrist_Gyroscope_X)    6.296  8.727  1.665  0.1034    
s(wrist_Gyroscope_Y)    3.662  4.528  2.729  0.0194 *  
s(wrist_Gyroscope_Z)    1.901  2.330  3.353  0.0326 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.538   Deviance explained = 58.7%
-REML score = 654.56  Scale est. = 1.0962    n = 430 (44 x 10)

model_accgyrx_sensors <- pffr(Y ~
                 +s( trunk_Accelerometer_X )
                 +s( upperarm_Accelerometer_X )
                 +s( wrist_Accelerometer_X )
                 +s( trunk_Gyroscope_X )
                 +s( upperarm_Gyroscope_X )
                 +s( wrist_Gyroscope_X ),
                 data=input_data_sensors, method = "REML", algorithm = "gam")

summary(model_accgyrx_sensors)

Family: gaussian 
Link function: identity 

Formula:
Y ~ +s(trunk_Accelerometer_X) + s(upperarm_Accelerometer_X) + 
    s(wrist_Accelerometer_X) + s(trunk_Gyroscope_X) + s(upperarm_Gyroscope_X) + 
    s(wrist_Gyroscope_X)

Constant coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.96009    0.04843   40.48   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Smooth terms & functional coefficients:
                               edf Ref.df      F  p-value    
Intercept(yindex)            7.116  9.000 45.040  < 2e-16 ***
s(trunk_Accelerometer_X)     1.009  1.016  3.518 0.059861 .  
s(upperarm_Accelerometer_X)  5.334  7.459  1.970 0.057583 .  
s(wrist_Accelerometer_X)    11.315 15.510  5.253  < 2e-16 ***
s(trunk_Gyroscope_X)        12.008 15.879  2.706 0.000459 ***
s(upperarm_Gyroscope_X)      1.457  1.852  0.935 0.326275    
s(wrist_Gyroscope_X)        14.566 18.681  3.042 2.89e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.618   Deviance explained = 66.5%
-REML score = 632.74  Scale est. = 0.90698   n = 430 (44 x 10)

model_accgyry_sensors <- pffr(Y ~
                 +s( trunk_Accelerometer_Y )
                 +s( upperarm_Accelerometer_Y )
                 +s( wrist_Accelerometer_Y )
                 +s( trunk_Gyroscope_Y )
                 +s( upperarm_Gyroscope_Y )
                 +s( wrist_Gyroscope_Y ),
                 data=input_data_sensors, method = "REML", algorithm = "gam")

summary(model_accgyry_sensors)

Family: gaussian 
Link function: identity 

Formula:
Y ~ +s(trunk_Accelerometer_Y) + s(upperarm_Accelerometer_Y) + 
    s(wrist_Accelerometer_Y) + s(trunk_Gyroscope_Y) + s(upperarm_Gyroscope_Y) + 
    s(wrist_Gyroscope_Y)

Constant coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   1.9313     0.0521   37.07   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Smooth terms & functional coefficients:
                              edf Ref.df      F  p-value    
Intercept(yindex)           6.579  9.000 33.862  < 2e-16 ***
s(trunk_Accelerometer_Y)    4.774  5.720  3.292 0.005144 ** 
s(upperarm_Accelerometer_Y) 4.618  5.532  1.543 0.313813    
s(wrist_Accelerometer_Y)    3.741  4.948  4.836 0.000268 ***
s(trunk_Gyroscope_Y)        3.329  4.857  1.420 0.221610    
s(upperarm_Gyroscope_Y)     4.447  5.616  1.862 0.083106 .  
s(wrist_Gyroscope_Y)        3.330  4.189  2.190 0.064472 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.511   Deviance explained = 54.6%
-REML score = 663.43  Scale est. = 1.1607    n = 430 (44 x 10)

model_accgyrz_sensors <- pffr(Y ~
                 +s( trunk_Accelerometer_Z )
                 +s( upperarm_Accelerometer_Z )
                 +s( wrist_Accelerometer_Z )
                 +s( trunk_Gyroscope_Z )
                 +s( upperarm_Gyroscope_Z )
                 +s( wrist_Gyroscope_Z ),
                 data=input_data_sensors, method = "REML", algorithm = "gam")

summary(model_accgyrz_sensors)

Family: gaussian 
Link function: identity 

Formula:
Y ~ +s(trunk_Accelerometer_Z) + s(upperarm_Accelerometer_Z) + 
    s(wrist_Accelerometer_Z) + s(trunk_Gyroscope_Z) + s(upperarm_Gyroscope_Z) + 
    s(wrist_Gyroscope_Z)

Constant coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.92772    0.04986   38.66   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Smooth terms & functional coefficients:
                               edf Ref.df      F  p-value    
Intercept(yindex)            7.024  9.000 40.527  < 2e-16 ***
s(trunk_Accelerometer_Z)     2.462  3.022  6.252 0.000363 ***
s(upperarm_Accelerometer_Z)  6.357  6.847 11.123  < 2e-16 ***
s(wrist_Accelerometer_Z)     1.000  1.001  5.156 0.023684 *  
s(trunk_Gyroscope_Z)         1.400  1.693  3.626 0.022934 *  
s(upperarm_Gyroscope_Z)     12.026 16.296  1.876 0.021428 *  
s(wrist_Gyroscope_Z)         1.001  1.001  9.151 0.002636 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.564   Deviance explained = 59.6%
-REML score = 643.27  Scale est. = 1.0351    n = 430 (44 x 10)

3.6 Downsampling the sensor signals

# Define the rate for down-sampling
sampling_rate = 10

# Define the standard length of a full 45 minute session (45min*60sec*100Hz = 27000)
standard_length <- 27000

decimate_stdz = function(df, sampling_rate , standard_length){
  
  # Define decimate function
  decimate_func = function(x) {signal::decimate(x, sampling_rate)}
  
  # Apply decimate to sensor signals columns and store in separate variables
  df2 <- lapply(df[,3:11], decimate_func)
  df2 = data.frame(
    setNames(df2, names(df)[3:11]))
  
  df3 = lapply(df2, function(x) {
    if (length(x) >= standard_length) {
      x <- x[1:standard_length]
    } else if (length(x) < standard_length) {
      x <- c(x, rep(0, standard_length - length(x)))
    }
    return(x)
  })
  
  # Downsample timestamp and id columns separately
  Timestamp_downsampled <- df$Timestamp[seq(1, nrow(df), sampling_rate)]
  sub_con_downsampled <- df$sub_con[seq(1, nrow(df), sampling_rate)]
  
  # Standardize Timestamp and sub_con
  if (nrow(df2) >= standard_length) {
    Timestamp <- Timestamp_downsampled[1:standard_length]
    sub_con <- sub_con_downsampled[1:standard_length]
  } else if (nrow(df2) < standard_length) {
    n_to_add <- standard_length - nrow(df2)
    Timestamp <- c(Timestamp_downsampled, Timestamp_downsampled[nrow(df2)] + (1/sampling_rate)*(1:n_to_add))
    sub_con <- c(sub_con_downsampled, rep(sub_con_downsampled[nrow(df2)], standard_length - nrow(df2)))
  }
  
  # Combine the results into a new data frame
  df3 <- data.frame(
    Timestamp = Timestamp,
    sub_con = sub_con,
    setNames(df3, names(df)[3:11]) # Ensure column names are preserved
  )
  
  list = list(decimate_std_df = df3
              )
  return(list)
}

trunk_total_df_dcmt <- trunk_total_df |> 
  dplyr::group_by(sub_con) |> 
  dplyr::do(decimate_stdz(. , standard_length = standard_length, sampling_rate = sampling_rate)$decimate_std_df) |>
  dplyr::select(!contains('Magnetometer')) |>
  dplyr::ungroup()

wrist_total_df_dcmt <- wrist_total_df |> 
  dplyr::group_by(sub_con) |> 
  dplyr::do(decimate_stdz(. , standard_length = standard_length, sampling_rate = sampling_rate)$decimate_std_df) |>
  dplyr::select(!contains('Magnetometer')) |>
  dplyr::ungroup()
  
upperarm_total_df_dcmt <- upperarm_total_df |> 
  dplyr::group_by(sub_con) |> 
  dplyr::do(decimate_stdz(. , standard_length = standard_length, sampling_rate = sampling_rate)$decimate_std_df) |>
  dplyr::select(!contains('Magnetometer')) |>
  dplyr::ungroup()

3.7 FPCA on all 30-second windows of the filtered (not downsampled) signals for the 5-minute intervals

3.7.1 Input preparation

# Define the pve (portion of variance explained)
pve = 0.99


# The function to prepare the dataframes for the FPCA function and apply the butterworth filter

df_prep_btrwrth = function(df, threshold=15, pve){
  df |> 
    dplyr::select(-contains('Magnetometer')) |>
    dplyr::mutate(time= as.POSIXct(df[,1], origin="1970-01-01")) |>
    dplyr::mutate_at(3:8,.funs=butter_filter_data) |>
    dplyr::group_by(sub_con) |>
    dplyr::mutate(time_diff_secs= as.numeric(difftime(time, time[1], units = "secs"))) |>
    dplyr::mutate(time_diff_mins= round(as.numeric(difftime(time, time[1], units = "mins")), digits = 0)) |>
    dplyr::mutate(max_time_diff_secs = max(time_diff_secs, na.rm = TRUE)) |>
    dplyr::mutate(interval = round(time_diff_secs / 300, digits = 0) * 300) |>
    dplyr::filter(abs(time_diff_secs - interval) <= threshold | (max_time_diff_secs-time_diff_secs)<=(2*threshold), time_diff_mins<=45)|>
    dplyr::mutate(interval = round(time_diff_mins / 5, digits = 0), index = seq(1, length(time))) |>
    dplyr::select(-max_time_diff_secs) |>
    dplyr::relocate(c(time,time_diff_secs,time_diff_mins,interval, index), .after = Timestamp)
  }
  
# The function to extract the principle components that explain at least up to pve that is specified

get_pcs_30s = function(df, pve) {
  set.seed(2022)
    col_names_adj = colnames(y_T_val)[colnames(y_T_val) %in% unique(df$sub_con)]
    final_result = matrix(1:length(col_names_adj))
    colnames(final_result) = 'sub_con'
    for (i in 1:(dim(df)[2]-7)){
      col = df[,c(6,7,(i+7))]
      temp_fpc_wide = pivot_wider(col, values_from = colnames(col)[3], names_from = index )
      temp_fpc_wide = as.data.frame(temp_fpc_wide)
      row.names(temp_fpc_wide) = temp_fpc_wide$sub_con
      temp_fpc_wide = temp_fpc_wide[,-1]
      temp_fpca_pace_model = refund::fpca.face(as.matrix(temp_fpc_wide), pve = pve)
      result = as.matrix(temp_fpca_pace_model[["scores"]])
      result <- as.matrix(result[order(match(rownames(result), col_names_adj)),])
      colnames(result) = paste0(colnames(col)[3],'_','pc',c(1:dim(result)[2]))
      final_result = cbind(final_result,  result)
    }
    final_result = as.data.frame(apply(final_result[,-1],2, as.numeric))
    rownames(final_result) = col_names_adj
    return(final_result)
  }

# Preparing the FPCA inputs for the model

trunk_total_df_pcs_30sec = get_pcs_30s(df_prep_btrwrth(trunk_total_df), pve)
upperarm_total_df_pcs_30sec = get_pcs_30s(df_prep_btrwrth(upperarm_total_df), pve)
wrist_total_df_pcs_30sec = get_pcs_30s(df_prep_btrwrth(wrist_total_df), pve)

# Input to the pffr model

input_data_fpca_30sec = input_data_sensors
input_data_fpca_30sec[10:36] <- NULL


for (i in 1:6){
  column = colnames(trunk_total_df)[(i+2)]
  df_temp <- trunk_total_df_pcs_30sec |> dplyr::select(contains(column))
  df_temp = as.matrix(df_temp)
  input_data_fpca_30sec[[(9+i)]] <-  input_data_fpca_30sec[[column]] <- df_temp
}

for (i in 1:6){
  column = colnames(wrist_total_df)[(i+2)]
  df_temp <- wrist_total_df_pcs_30sec |> dplyr::select(contains(column))
  df_temp <- rbind(df_temp[1:2,],c(rep(1e-6,(dim(df_temp)[2]))),df_temp[3:43,])
  df_temp = as.matrix(df_temp)
  rownames(df_temp) = colnames(y_T_val)
  input_data_fpca_30sec[[(15+i)]] <-  input_data_fpca_30sec[[column]] <- df_temp
}

for (i in 1:6){
  column = colnames(upperarm_total_df)[(i+2)]
  df_temp <- upperarm_total_df_pcs_30sec |> dplyr::select(contains(column))
  df_temp = as.matrix(df_temp)
  input_data_fpca_30sec[[(21+i)]] <-  input_data_fpca_30sec[[column]] <- df_temp
}

3.7.2 Modeling

set.seed(2022)

model_acc_pcs_30sec <- pffr(Y ~
                 c(xload)  +
                 c(xpace)  +  
                 c(xgender)+ 
                 c(xage)+
                 c(xheight)+ 
                 c(xweight)+ 
                 c(xhip)+ 
                 c(xwaist)+
                ff( trunk_Accelerometer_X )+
                ff( trunk_Accelerometer_Y )+
                ff( trunk_Accelerometer_Z )+
                ff( upperarm_Accelerometer_X )+
                ff( upperarm_Accelerometer_Y )+
                ff( upperarm_Accelerometer_Z )+
                ff( wrist_Accelerometer_X )+
                ff( wrist_Accelerometer_Y )+
                ff( wrist_Accelerometer_Z ),
                 data=input_data_fpca_30sec , method = "REML", algorithm = "gam")

summary(model_acc_pcs_30sec)

Family: gaussian 
Link function: identity 

Formula:
Y ~ c(xload) + c(xpace) + c(xgender) + c(xage) + c(xheight) + 
    c(xweight) + c(xhip) + c(xwaist) + ff(trunk_Accelerometer_X) + 
    ff(trunk_Accelerometer_Y) + ff(trunk_Accelerometer_Z) + ff(upperarm_Accelerometer_X) + 
    ff(upperarm_Accelerometer_Y) + ff(upperarm_Accelerometer_Z) + 
    ff(wrist_Accelerometer_X) + ff(wrist_Accelerometer_Y) + ff(wrist_Accelerometer_Z)

Constant coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -24.762956   2.952456  -8.387 1.86e-15 ***
xload         1.909958   0.216952   8.804  < 2e-16 ***
xpace         0.197673   0.024029   8.226 5.62e-15 ***
xgenderM     -3.121010   0.734748  -4.248 2.87e-05 ***
xage         -0.069726   0.039632  -1.759 0.079521 .  
xheight       0.072898   0.020188   3.611 0.000356 ***
xweight      -0.009612   0.006154  -1.562 0.119307    
xhip          0.073791   0.015452   4.776 2.79e-06 ***
xwaist        0.054662   0.020165   2.711 0.007093 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Smooth terms & functional coefficients:
                                edf Ref.df       F p-value    
Intercept(yindex)             8.643  9.000 388.332  <2e-16 ***
ff(trunk_Accelerometer_X)     8.351  9.861   7.316  <2e-16 ***
ff(trunk_Accelerometer_Y)    16.038 18.273   7.250  <2e-16 ***
ff(trunk_Accelerometer_Z)    11.781 13.664  22.591  <2e-16 ***
ff(upperarm_Accelerometer_X) 13.955 15.914  12.135  <2e-16 ***
ff(upperarm_Accelerometer_Y) 13.555 15.594  10.659  <2e-16 ***
ff(upperarm_Accelerometer_Z) 17.119 19.169  16.308  <2e-16 ***
ff(wrist_Accelerometer_X)    11.627 13.433   8.830  <2e-16 ***
ff(wrist_Accelerometer_Y)    11.368 13.438  19.768  <2e-16 ***
ff(wrist_Accelerometer_Z)    12.977 14.991   7.468  <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.951   Deviance explained = 96.6%
-REML score = 416.55  Scale est. = 0.11547   n = 440 (44 x 10)

model_acc2_pcs_30sec <- pffr(Y ~
                ff( trunk_Accelerometer_X )+
                ff( trunk_Accelerometer_Y )+
                ff( trunk_Accelerometer_Z )+
                ff( upperarm_Accelerometer_X )+
                ff( upperarm_Accelerometer_Y )+
                ff( upperarm_Accelerometer_Z )+
                ff( wrist_Accelerometer_X )+
                ff( wrist_Accelerometer_Y )+
                ff( wrist_Accelerometer_Z ),
                 data=input_data_fpca_30sec , method = "REML", algorithm = "gam")

summary(model_acc2_pcs_30sec)

Family: gaussian 
Link function: identity 

Formula:
Y ~ ff(trunk_Accelerometer_X) + ff(trunk_Accelerometer_Y) + ff(trunk_Accelerometer_Z) + 
    ff(upperarm_Accelerometer_X) + ff(upperarm_Accelerometer_Y) + 
    ff(upperarm_Accelerometer_Z) + ff(wrist_Accelerometer_X) + 
    ff(wrist_Accelerometer_Y) + ff(wrist_Accelerometer_Z)

Constant coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.87322    0.02265    82.7   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Smooth terms & functional coefficients:
                                edf Ref.df       F p-value    
Intercept(yindex)             8.364  9.000 208.222  <2e-16 ***
ff(trunk_Accelerometer_X)     8.071  9.963   7.999  <2e-16 ***
ff(trunk_Accelerometer_Y)     4.448  4.749  48.304  <2e-16 ***
ff(trunk_Accelerometer_Z)    10.940 12.682  17.724  <2e-16 ***
ff(upperarm_Accelerometer_X) 11.149 12.730  33.860  <2e-16 ***
ff(upperarm_Accelerometer_Y) 12.945 15.078  13.094  <2e-16 ***
ff(upperarm_Accelerometer_Z) 13.850 15.695  69.436  <2e-16 ***
ff(wrist_Accelerometer_X)     9.023 10.599  33.496  <2e-16 ***
ff(wrist_Accelerometer_Y)     2.737  3.082  87.695  <2e-16 ***
ff(wrist_Accelerometer_Z)    10.217 11.807  35.507  <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.908   Deviance explained = 92.7%
-REML score = 467.36  Scale est. = 0.21771   n = 440 (44 x 10)

model_accgyrx_pcs_30sec <- pffr(Y ~
                ff( trunk_Accelerometer_X )+
                ff( upperarm_Accelerometer_X )+
                ff( wrist_Accelerometer_X )+
                ff( trunk_Gyroscope_X )+
                ff( upperarm_Gyroscope_X )+
                ff( wrist_Gyroscope_X ),
                 data=input_data_fpca_30sec , method = "REML", algorithm = "gam")

summary(model_accgyrx_pcs_30sec)

Family: gaussian 
Link function: identity 

Formula:
Y ~ ff(trunk_Accelerometer_X) + ff(upperarm_Accelerometer_X) + 
    ff(wrist_Accelerometer_X) + ff(trunk_Gyroscope_X) + ff(upperarm_Gyroscope_X) + 
    ff(wrist_Gyroscope_X)

Constant coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.90348    0.03939   48.33   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Smooth terms & functional coefficients:
                               edf Ref.df      F  p-value    
Intercept(yindex)            7.422  9.000 67.525  < 2e-16 ***
ff(trunk_Accelerometer_X)    1.960  2.333  2.288   0.0991 .  
ff(upperarm_Accelerometer_X) 8.084  9.725  6.115  < 2e-16 ***
ff(wrist_Accelerometer_X)    9.094 10.714  9.759  < 2e-16 ***
ff(trunk_Gyroscope_X)        6.682  7.708 19.512  < 2e-16 ***
ff(upperarm_Gyroscope_X)     4.747  4.933 13.741  < 2e-16 ***
ff(wrist_Gyroscope_X)        6.813  8.569  4.320 3.76e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.714   Deviance explained = 74.3%
-REML score = 626.71  Scale est. = 0.67343   n = 440 (44 x 10)

model_accgyry_pcs_30sec <- pffr(Y ~
                ff( trunk_Accelerometer_Y )+
                ff( upperarm_Accelerometer_Y )+
                ff( wrist_Accelerometer_Y )+
                ff( trunk_Gyroscope_Y )+
                ff( upperarm_Gyroscope_Y )+
                ff( wrist_Gyroscope_Y ),
                 data=input_data_fpca_30sec , method = "REML", algorithm = "gam")

summary(model_accgyry_pcs_30sec)

Family: gaussian 
Link function: identity 

Formula:
Y ~ ff(trunk_Accelerometer_Y) + ff(upperarm_Accelerometer_Y) + 
    ff(wrist_Accelerometer_Y) + ff(trunk_Gyroscope_Y) + ff(upperarm_Gyroscope_Y) + 
    ff(wrist_Gyroscope_Y)

Constant coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.91903    0.03762   51.01   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Smooth terms & functional coefficients:
                               edf Ref.df      F  p-value    
Intercept(yindex)            7.522  9.000 74.333  < 2e-16 ***
ff(trunk_Accelerometer_Y)    9.021 10.542 21.694  < 2e-16 ***
ff(upperarm_Accelerometer_Y) 6.665  8.301 11.992  < 2e-16 ***
ff(wrist_Accelerometer_Y)    8.333 10.014  8.909  < 2e-16 ***
ff(trunk_Gyroscope_Y)        7.950  9.441 10.091  < 2e-16 ***
ff(upperarm_Gyroscope_Y)     7.020  8.309  5.435 1.47e-06 ***
ff(wrist_Gyroscope_Y)        3.312  3.649 32.716  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.738   Deviance explained = 76.8%
-REML score = 616.78  Scale est. = 0.61753   n = 440 (44 x 10)

model_accgyrz_pcs_30sec <- pffr(Y ~
                ff( trunk_Accelerometer_Z )+
                ff( upperarm_Accelerometer_Z )+
                ff( wrist_Accelerometer_Z )+
                ff( trunk_Gyroscope_Z )+
                ff( upperarm_Gyroscope_Z )+
                ff( wrist_Gyroscope_Z ),
                 data=input_data_fpca_30sec , method = "REML", algorithm = "gam")

summary(model_accgyrz_pcs_30sec)

Family: gaussian 
Link function: identity 

Formula:
Y ~ ff(trunk_Accelerometer_Z) + ff(upperarm_Accelerometer_Z) + 
    ff(wrist_Accelerometer_Z) + ff(trunk_Gyroscope_Z) + ff(upperarm_Gyroscope_Z) + 
    ff(wrist_Gyroscope_Z)

Constant coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.88403    0.03583   52.58   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Smooth terms & functional coefficients:
                                edf Ref.df      F  p-value    
Intercept(yindex)             7.618  9.000 81.081  < 2e-16 ***
ff(trunk_Accelerometer_Z)     6.550  8.291  5.812 7.24e-07 ***
ff(upperarm_Accelerometer_Z) 10.224 12.073  9.905  < 2e-16 ***
ff(wrist_Accelerometer_Z)     3.992  4.278  6.894 2.81e-05 ***
ff(trunk_Gyroscope_Z)         7.581  9.019  8.359  < 2e-16 ***
ff(upperarm_Gyroscope_Z)      1.524  1.820 10.434 6.33e-05 ***
ff(wrist_Gyroscope_Z)        10.161 12.306  7.812  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.763   Deviance explained = 78.9%
-REML score = 577.85  Scale est. = 0.55906   n = 440 (44 x 10)

model_gyr_pcs_30sec <- pffr(Y ~
                 c(xload)  +
                 c(xpace)  +  
                 c(xgender)+ 
                 c(xage)+
                 c(xheight)+ 
                 c(xweight)+ 
                 c(xhip)+ 
                 c(xwaist)+
                 ff( trunk_Gyroscope_X )+
                 ff( trunk_Gyroscope_Y )+
                 ff( trunk_Gyroscope_Z )+
                 ff( upperarm_Gyroscope_X )+
                 ff( upperarm_Gyroscope_Y )+
                 ff( upperarm_Gyroscope_Z )+
                 ff( wrist_Gyroscope_X )+
                 ff( wrist_Gyroscope_Y )+
                 ff( wrist_Gyroscope_Z ),
                 data=input_data_fpca_30sec , method = "REML", algorithm = "gam")

summary(model_gyr_pcs_30sec)

Family: gaussian 
Link function: identity 

Formula:
Y ~ c(xload) + c(xpace) + c(xgender) + c(xage) + c(xheight) + 
    c(xweight) + c(xhip) + c(xwaist) + ff(trunk_Gyroscope_X) + 
    ff(trunk_Gyroscope_Y) + ff(trunk_Gyroscope_Z) + ff(upperarm_Gyroscope_X) + 
    ff(upperarm_Gyroscope_Y) + ff(upperarm_Gyroscope_Z) + ff(wrist_Gyroscope_X) + 
    ff(wrist_Gyroscope_Y) + ff(wrist_Gyroscope_Z)

Constant coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.8689657  3.1777068   1.218 0.224295    
xload       -0.1970625  0.1824696  -1.080 0.280964    
xpace       -0.0247752  0.0210811  -1.175 0.240773    
xgenderM    -0.1501674  0.5475435  -0.274 0.784064    
xage         0.0463797  0.0263344   1.761 0.079159 .  
xheight     -0.0086992  0.0205231  -0.424 0.671943    
xweight      0.0207541  0.0053317   3.893 0.000121 ***
xhip        -0.0253390  0.0069905  -3.625 0.000336 ***
xwaist       0.0006995  0.0088369   0.079 0.936961    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Smooth terms & functional coefficients:
                            edf Ref.df       F p-value    
Intercept(yindex)         8.626  9.000 357.243 < 2e-16 ***
ff(trunk_Gyroscope_X)    18.598 20.439  22.449 < 2e-16 ***
ff(trunk_Gyroscope_Y)     2.899  3.041  13.070 < 2e-16 ***
ff(trunk_Gyroscope_Z)    12.430 14.781   9.471 < 2e-16 ***
ff(upperarm_Gyroscope_X) 13.957 15.943  11.917 < 2e-16 ***
ff(upperarm_Gyroscope_Y) 12.714 15.050  10.505 < 2e-16 ***
ff(upperarm_Gyroscope_Z)  5.383  6.649   6.268 1.7e-06 ***
ff(wrist_Gyroscope_X)    12.812 14.756  10.001 < 2e-16 ***
ff(wrist_Gyroscope_Y)     6.342  7.744   2.032  0.0411 *  
ff(wrist_Gyroscope_Z)    16.047 18.291  31.796 < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.949   Deviance explained = 96.3%
-REML score = 420.41  Scale est. = 0.11912   n = 440 (44 x 10)

model_gyr2_pcs_30sec <- pffr(Y ~
                 ff( trunk_Gyroscope_X )+
                 ff( trunk_Gyroscope_Y )+
                 ff( trunk_Gyroscope_Z )+
                 ff( upperarm_Gyroscope_X )+
                 ff( upperarm_Gyroscope_Y )+
                 ff( upperarm_Gyroscope_Z )+
                 ff( wrist_Gyroscope_X )+
                 ff( wrist_Gyroscope_Y )+
                 ff( wrist_Gyroscope_Z ),
                 data=input_data_fpca_30sec , method = "REML", algorithm = "gam")

summary(model_gyr2_pcs_30sec)

Family: gaussian 
Link function: identity 

Formula:
Y ~ ff(trunk_Gyroscope_X) + ff(trunk_Gyroscope_Y) + ff(trunk_Gyroscope_Z) + 
    ff(upperarm_Gyroscope_X) + ff(upperarm_Gyroscope_Y) + ff(upperarm_Gyroscope_Z) + 
    ff(wrist_Gyroscope_X) + ff(wrist_Gyroscope_Y) + ff(wrist_Gyroscope_Z)

Constant coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.88076    0.01768   106.4   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Smooth terms & functional coefficients:
                            edf Ref.df       F  p-value    
Intercept(yindex)         8.601  9.000 335.431  < 2e-16 ***
ff(trunk_Gyroscope_X)    18.332 20.195  47.453  < 2e-16 ***
ff(trunk_Gyroscope_Y)     4.636  4.855  25.906  < 2e-16 ***
ff(trunk_Gyroscope_Z)    12.322 14.861  24.293  < 2e-16 ***
ff(upperarm_Gyroscope_X) 14.324 16.505  27.992  < 2e-16 ***
ff(upperarm_Gyroscope_Y) 12.911 15.334  25.105  < 2e-16 ***
ff(upperarm_Gyroscope_Z)  5.819  7.244   6.118 1.07e-06 ***
ff(wrist_Gyroscope_X)    12.768 14.613  25.653  < 2e-16 ***
ff(wrist_Gyroscope_Y)     3.390  4.386   1.699     0.14    
ff(wrist_Gyroscope_Z)    16.008 18.218  52.384  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.946   Deviance explained = 95.9%
-REML score = 410.13  Scale est. = 0.12787   n = 440 (44 x 10)

model_trunk_pcs_30sec <- pffr(Y ~
                 c(xload)  +
                 c(xpace)  +  
                 c(xgender)+ 
                 c(xage)+
                 c(xheight)+ 
                 c(xweight)+ 
                 c(xhip)+ 
                 c(xwaist)+
                 ff( trunk_Accelerometer_X )+
                 ff( trunk_Accelerometer_Y )+
                 ff( trunk_Accelerometer_Z )+
                 ff( trunk_Gyroscope_X )+
                 ff( trunk_Gyroscope_Y )+
                 ff( trunk_Gyroscope_Z ),
                 data=input_data_fpca_30sec , method = "REML", algorithm = "gam")

summary(model_trunk_pcs_30sec)

Family: gaussian 
Link function: identity 

Formula:
Y ~ c(xload) + c(xpace) + c(xgender) + c(xage) + c(xheight) + 
    c(xweight) + c(xhip) + c(xwaist) + ff(trunk_Accelerometer_X) + 
    ff(trunk_Accelerometer_Y) + ff(trunk_Accelerometer_Z) + ff(trunk_Gyroscope_X) + 
    ff(trunk_Gyroscope_Y) + ff(trunk_Gyroscope_Z)

Constant coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -4.229881   1.198135  -3.530 0.000466 ***
xload        0.938783   0.124954   7.513 4.21e-13 ***
xpace        0.136999   0.013251  10.339  < 2e-16 ***
xgenderM    -0.281628   0.269111  -1.047 0.295993    
xage         0.028977   0.020461   1.416 0.157545    
xheight     -0.029801   0.009528  -3.128 0.001898 ** 
xweight      0.004502   0.005417   0.831 0.406448    
xhip         0.080728   0.007745  10.423  < 2e-16 ***
xwaist      -0.017332   0.006846  -2.532 0.011752 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Smooth terms & functional coefficients:
                            edf Ref.df       F  p-value    
Intercept(yindex)         7.989  9.000 121.502  < 2e-16 ***
ff(trunk_Accelerometer_X) 5.791  7.378   4.000 0.000305 ***
ff(trunk_Accelerometer_Y) 7.485  9.387   4.505 9.88e-06 ***
ff(trunk_Accelerometer_Z) 9.783 11.800   9.926  < 2e-16 ***
ff(trunk_Gyroscope_X)     5.735  6.551  11.189  < 2e-16 ***
ff(trunk_Gyroscope_Y)     8.569 10.168   5.247 2.82e-07 ***
ff(trunk_Gyroscope_Z)     7.292  9.071   5.789  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.839   Deviance explained = 86.2%
-REML score = 522.11  Scale est. = 0.37836   n = 440 (44 x 10)

model_wrist_pcs_30sec <- pffr(Y ~
                 c(xload)  +
                 c(xpace)  +  
                 c(xgender)+ 
                 c(xage)+
                 c(xheight)+ 
                 c(xweight)+ 
                 c(xhip)+ 
                 c(xwaist)+
                 ff( wrist_Accelerometer_X )+
                 ff( wrist_Accelerometer_Y )+
                 ff( wrist_Accelerometer_Z )+
                 ff( wrist_Gyroscope_X )+
                 ff( wrist_Gyroscope_Y )+
                 ff( wrist_Gyroscope_Z ),
                 data=input_data_fpca_30sec , method = "REML", algorithm = "gam")

summary(model_wrist_pcs_30sec)

Family: gaussian 
Link function: identity 

Formula:
Y ~ c(xload) + c(xpace) + c(xgender) + c(xage) + c(xheight) + 
    c(xweight) + c(xhip) + c(xwaist) + ff(wrist_Accelerometer_X) + 
    ff(wrist_Accelerometer_Y) + ff(wrist_Accelerometer_Z) + ff(wrist_Gyroscope_X) + 
    ff(wrist_Gyroscope_Y) + ff(wrist_Gyroscope_Z)

Constant coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -13.957035   1.118037 -12.484  < 2e-16 ***
xload         1.509295   0.159079   9.488  < 2e-16 ***
xpace         0.308940   0.027554  11.212  < 2e-16 ***
xgenderM     -0.082014   0.280899  -0.292 0.770475    
xage         -0.122113   0.029135  -4.191 3.48e-05 ***
xheight       0.027667   0.007725   3.581 0.000388 ***
xweight      -0.076963   0.007632 -10.084  < 2e-16 ***
xhip          0.128497   0.008477  15.159  < 2e-16 ***
xwaist        0.005042   0.005705   0.884 0.377348    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Smooth terms & functional coefficients:
                             edf Ref.df       F p-value    
Intercept(yindex)          8.269  9.000 178.400  <2e-16 ***
ff(wrist_Accelerometer_X)  7.422  8.916  15.421  <2e-16 ***
ff(wrist_Accelerometer_Y)  8.105  9.839  18.287  <2e-16 ***
ff(wrist_Accelerometer_Z)  4.502  4.734  17.426  <2e-16 ***
ff(wrist_Gyroscope_X)     11.140 12.956   9.395  <2e-16 ***
ff(wrist_Gyroscope_Y)     10.835 12.511  19.904  <2e-16 ***
ff(wrist_Gyroscope_Z)     13.953 15.847   9.016  <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.891   Deviance explained = 90.9%
-REML score = 493.75  Scale est. = 0.25794   n = 440 (44 x 10)

model_upperarm_pcs_30sec <- pffr(Y ~
                 c(xload)  +
                 c(xpace)  +  
                 c(xgender)+ 
                 c(xage)+
                 c(xheight)+ 
                 c(xweight)+ 
                 c(xhip)+ 
                 c(xwaist)+
                 ff( upperarm_Accelerometer_X )+
                 ff( upperarm_Accelerometer_Y )+
                 ff( upperarm_Accelerometer_Z )+
                 ff( upperarm_Gyroscope_X )+
                 ff( upperarm_Gyroscope_Y )+
                 ff( upperarm_Gyroscope_Z ),
                 data=input_data_fpca_30sec , method = "REML", algorithm = "gam")

summary(model_upperarm_pcs_30sec)

Family: gaussian 
Link function: identity 

Formula:
Y ~ c(xload) + c(xpace) + c(xgender) + c(xage) + c(xheight) + 
    c(xweight) + c(xhip) + c(xwaist) + ff(upperarm_Accelerometer_X) + 
    ff(upperarm_Accelerometer_Y) + ff(upperarm_Accelerometer_Z) + 
    ff(upperarm_Gyroscope_X) + ff(upperarm_Gyroscope_Y) + ff(upperarm_Gyroscope_Z)

Constant coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.866884   1.969823   0.440 0.660126    
xload        0.326676   0.141406   2.310 0.021408 *  
xpace        0.161836   0.013705  11.808  < 2e-16 ***
xgenderM     1.609993   0.476617   3.378 0.000805 ***
xage         0.099186   0.022298   4.448 1.14e-05 ***
xheight     -0.039214   0.012838  -3.054 0.002412 ** 
xweight     -0.011814   0.006549  -1.804 0.072019 .  
xhip         0.048308   0.008406   5.747 1.85e-08 ***
xwaist      -0.028376   0.011007  -2.578 0.010306 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Smooth terms & functional coefficients:
                                edf Ref.df       F  p-value    
Intercept(yindex)             7.959  9.000 116.753  < 2e-16 ***
ff(upperarm_Accelerometer_X)  9.665 11.533  10.589  < 2e-16 ***
ff(upperarm_Accelerometer_Y)  1.980  2.371   1.848    0.126    
ff(upperarm_Accelerometer_Z) 10.493 12.441  12.368  < 2e-16 ***
ff(upperarm_Gyroscope_X)      3.208  3.507  11.875 1.96e-07 ***
ff(upperarm_Gyroscope_Y)      5.714  6.958   7.378  < 2e-16 ***
ff(upperarm_Gyroscope_Z)      8.575 10.357   9.339  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.835   Deviance explained = 85.6%
-REML score = 526.12  Scale est. = 0.38976   n = 440 (44 x 10)

3.8 Summary of Results For the Functional Models

models <- list(model_1, model_2, model_mad, model_max, model_mean, model_med, model_min, model_Perc25, model_Perc75, model_rms,model_std, model_trunk_sensors, model_upperarm_sensors, model_wrist_sensors,model_acc_sensors,model_gyr_sensors,model_acc2_sensors,model_gyr2_sensors, model_accgyrx_sensors,model_accgyry_sensors,model_accgyrz_sensors,model_acc_pcs_30sec,model_gyr_pcs_30sec,model_trunk_pcs_30sec,model_wrist_pcs_30sec,model_upperarm_pcs_30sec,model_acc2_pcs_30sec,model_gyr2_pcs_30sec, model_accgyrx_pcs_30sec,model_accgyry_pcs_30sec, model_accgyrz_pcs_30sec)


#functions
rsq_adj = function(x) {return(summary(x)$r.sq)}
rsq = function(x) {return(summary(x)$dev)}
aic = function(x) {return(x$aic)}

summary_df <- data.frame(
  Model = c('model_task', 'model_task&prtcpntsChar', 'model_mad', 'model_max', 'model_mean','model_med', 'model_min', 'model_Perc25', 'model_Perc75', 'model_rms','model_std','model_trunk_sensors', 'model_upperarm_sensors', 'model_wrist_sensors','model_acc_sensors','model_gyr_sensors', 'model_acc2_sensors','model_gyr2_sensors', 'model_accgyrx_sensors','model_accgyry_sensors','model_accgyrz_sensors',
            'model_acc_pcs_30sec','model_gyr_pcs_30sec','model_trunk_pcs_30sec','model_wrist_pcs_30sec','model_upperarm_pcs_30sec', 'model_acc2_pcs_30sec','model_gyr2_pcs_30sec', 'model_accgyrx_pcs_30sec','model_accgyry_pcs_30sec', 'model_accgyrz_pcs_30sec'), 
  R_squared_adjusted = sapply(models, rsq_adj),
  R_squared = sapply(models, rsq),
  AIC = sapply(models, aic)
) |> dplyr::arrange(Model)
  #dplyr::arrange(desc(R_squared_adjusted))


# printing the summary of results

DT::datatable(
  data = summary_df, rownames = F, 
  extensions = c('Buttons','FixedColumns'),
  options = list(
    dom = 'Bfrtip',
    buttons = c('csv', 'excel', 'pdf', 'print'),
    paging = TRUE, searching = TRUE, info = FALSE,
    sort = TRUE, scrollX = TRUE, fixedColumns = list(leftColumns = 1)
  )
) |> 
  DT::formatRound(
    columns = c('R_squared_adjusted','R_squared', 'AIC'), 
    digits = 3)

4 Benchmark Models

We employed linear regression to (i) investigate the relationships in question from a cross-sectional perspective, similar to previous literature, and (ii) use these models as benchmarks to assess whether functional regression models are better at explaining the variability in upper extremity fatigue development.

4.1 Input

In the code chunk below, we aggregated the data across all time markers (i.e., over 45 minutes at time 0, 5, 10, …, 45) for each statistical feature as well as the response variable TRPE.

## Aggregate the data
agg_df = input_data_nested |>
    dplyr::mutate(across(-c(1:10), \(x) purrr::map_dbl(x, mean)))

colnames(agg_df)[c(11:65)] = paste0("mean_", colnames(agg_df)[c(11:65)])

DT::datatable(data = agg_df, rownames = F, class = list("display", "compact"), extensions = "FixedColumns",
    options = list(dom = "Bfrtip", paging = TRUE, info = FALSE, sort = TRUE, scrollX = TRUE,
        fixedColumns = list(leftColumns = 1))) |>
    formatRound(digits = 2, columns = c(11:65))

4.2 Modeling

4.2.1 With only task specifications

lm_task = lm(mean_TRPE ~ xpace + xload, agg_df)
summary(lm_task)

Call:
lm(formula = mean_TRPE ~ xpace + xload, data = agg_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.2665 -0.5966 -0.1763  0.6287  2.1384 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept) -0.68991    1.25461  -0.550   0.5860  
xpace        0.63447    0.40664   1.560   0.1280  
xload        0.08898    0.04368   2.037   0.0495 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9266 on 34 degrees of freedom
Multiple R-squared:  0.1159,    Adjusted R-squared:  0.06385 
F-statistic: 2.228 on 2 and 34 DF,  p-value: 0.1233

4.2.2 With only individual characteristics

lm_ind = lm(mean_TRPE ~ xgender + xage + xheight + xweight + xwaist + xhip, agg_df) 
summary(lm_ind)

Call:
lm(formula = mean_TRPE ~ xgender + xage + xheight + xweight + 
    xwaist + xhip, data = agg_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.1096 -0.4641 -0.1715  0.3495  1.7632 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  1.32994    5.55712   0.239  0.81248   
xgenderF    -0.90293    1.17846  -0.766  0.44955   
xage         0.06579    0.05687   1.157  0.25645   
xheight     -0.02923    0.03407  -0.858  0.39773   
xweight      0.01071    0.01749   0.613  0.54475   
xwaist      -0.05147    0.03437  -1.498  0.14462   
xhip         0.07600    0.02662   2.855  0.00774 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7559 on 30 degrees of freedom
Multiple R-squared:  0.4809,    Adjusted R-squared:  0.3771 
F-statistic: 4.632 on 6 and 30 DF,  p-value: 0.001927

4.2.3 Task specifications + individual characteristics

lm_task_ind = lm(mean_TRPE ~ xpace + xload + xgender + xage + xheight + xweight + xwaist + xhip, agg_df)

summary(lm_task_ind)

Call:
lm(formula = mean_TRPE ~ xpace + xload + xgender + xage + xheight + 
    xweight + xwaist + xhip, data = agg_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.8882 -0.3827 -0.1971  0.2437  1.8510 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  0.66238    4.86546   0.136  0.89268   
xpace        0.78982    0.29047   2.719  0.01111 * 
xload        0.10144    0.03163   3.207  0.00335 **
xgenderF    -1.37499    1.03423  -1.329  0.19442   
xage         0.08676    0.04984   1.741  0.09272 . 
xheight     -0.04439    0.03000  -1.480  0.15014   
xweight      0.01499    0.01527   0.982  0.33457   
xwaist      -0.06287    0.03008  -2.090  0.04578 * 
xhip         0.08220    0.02326   3.535  0.00144 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.657 on 28 degrees of freedom
Multiple R-squared:  0.6339,    Adjusted R-squared:  0.5293 
F-statistic: 6.061 on 8 and 28 DF,  p-value: 0.0001512

4.2.4 Adding aggregated features

lm_features = lm(mean_TRPE ~ xpace + xload + xgender + xage + xheight + xweight + xwaist + xhip + 
                   mean_Wrist_Acc_Mean + 
                   mean_Wrist_Jerk_Mean + 
                   mean_Trunk_Acc_Mean + 
                   mean_Trunk_Jerk_Mean +
                   mean_Upperarm_Acc_Mean + 
                   mean_Upperarm_Jerk_Mean + 
                   mean_Wrist_Acc_Std + 
                   mean_Wrist_Jerk_Std +
                   mean_Trunk_Acc_Std + 
                   mean_Trunk_Jerk_Std +
                   mean_Upperarm_Acc_Std + 
                   mean_Upperarm_Jerk_Std,
    agg_df)

summary(lm_features)

Call:
lm(formula = mean_TRPE ~ xpace + xload + xgender + xage + xheight + 
    xweight + xwaist + xhip + mean_Wrist_Acc_Mean + mean_Wrist_Jerk_Mean + 
    mean_Trunk_Acc_Mean + mean_Trunk_Jerk_Mean + mean_Upperarm_Acc_Mean + 
    mean_Upperarm_Jerk_Mean + mean_Wrist_Acc_Std + mean_Wrist_Jerk_Std + 
    mean_Trunk_Acc_Std + mean_Trunk_Jerk_Std + mean_Upperarm_Acc_Std + 
    mean_Upperarm_Jerk_Std, data = agg_df)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.61298 -0.21572 -0.01331  0.18417  0.84100 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)   
(Intercept)             -31.32659   25.08549  -1.249  0.22970   
xpace                     1.24089    0.31831   3.898  0.00128 **
xload                    -0.05246    0.07791  -0.673  0.51031   
xgenderF                 -2.19634    1.31736  -1.667  0.11492   
xage                      0.14275    0.05888   2.424  0.02754 * 
xheight                  -0.07274    0.03963  -1.836  0.08506 . 
xweight                   0.03816    0.02392   1.595  0.13022   
xwaist                   -0.11748    0.04212  -2.789  0.01313 * 
xhip                      0.09232    0.02543   3.631  0.00225 **
mean_Wrist_Acc_Mean      40.33907   23.28770   1.732  0.10247   
mean_Wrist_Jerk_Mean     -1.06247    0.69401  -1.531  0.14532   
mean_Trunk_Acc_Mean      14.79602   10.63697   1.391  0.18327   
mean_Trunk_Jerk_Mean      5.28882    3.27783   1.614  0.12618   
mean_Upperarm_Acc_Mean  -18.02929   11.65061  -1.547  0.14129   
mean_Upperarm_Jerk_Mean   0.86136    1.17385   0.734  0.47370   
mean_Wrist_Acc_Std       14.11491    6.84803   2.061  0.05593 . 
mean_Wrist_Jerk_Std      -0.04953    0.50185  -0.099  0.92261   
mean_Trunk_Acc_Std      -83.71895   71.09408  -1.178  0.25618   
mean_Trunk_Jerk_Std      -1.24167    1.37436  -0.903  0.37969   
mean_Upperarm_Acc_Std    11.09451   20.17185   0.550  0.58992   
mean_Upperarm_Jerk_Std   -1.19405    0.81594  -1.463  0.16272   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.5421 on 16 degrees of freedom
Multiple R-squared:  0.8576,    Adjusted R-squared:  0.6796 
F-statistic: 4.817 on 20 and 16 DF,  p-value: 0.001229

lm_features_1 = lm(mean_TRPE ~ xpace + xload + xgender + xage + xheight + xweight + xwaist + xhip + 
                     mean_Wrist_Acc_Mean + 
                     mean_Wrist_Jerk_Mean + 
                     mean_Trunk_Acc_Mean + 
                     mean_Trunk_Jerk_Mean +
                     mean_Upperarm_Acc_Mean + 
                     mean_Upperarm_Jerk_Mean,
    agg_df)

summary(lm_features_1)

Call:
lm(formula = mean_TRPE ~ xpace + xload + xgender + xage + xheight + 
    xweight + xwaist + xhip + mean_Wrist_Acc_Mean + mean_Wrist_Jerk_Mean + 
    mean_Trunk_Acc_Mean + mean_Trunk_Jerk_Mean + mean_Upperarm_Acc_Mean + 
    mean_Upperarm_Jerk_Mean, data = agg_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.7579 -0.3911 -0.2073  0.4172  1.4036 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)             -29.18426   24.16026  -1.208 0.239894    
xpace                     0.91704    0.32270   2.842 0.009490 ** 
xload                     0.07030    0.06040   1.164 0.256983    
xgenderF                 -2.89675    1.23490  -2.346 0.028411 *  
xage                      0.12043    0.05813   2.072 0.050217 .  
xheight                  -0.08942    0.03721  -2.403 0.025127 *  
xweight                   0.01868    0.01939   0.963 0.345848    
xwaist                   -0.09637    0.03472  -2.776 0.011021 *  
xhip                      0.11425    0.02695   4.239 0.000336 ***
mean_Wrist_Acc_Mean      29.69756   20.23604   1.468 0.156376    
mean_Wrist_Jerk_Mean     -0.07270    0.38224  -0.190 0.850894    
mean_Trunk_Acc_Mean      10.29152   11.05319   0.931 0.361913    
mean_Trunk_Jerk_Mean      0.04869    0.54274   0.090 0.929335    
mean_Upperarm_Acc_Mean   -3.79883   11.30471  -0.336 0.740025    
mean_Upperarm_Jerk_Mean  -0.09564    0.42406  -0.226 0.823636    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6539 on 22 degrees of freedom
Multiple R-squared:  0.7151,    Adjusted R-squared:  0.5338 
F-statistic: 3.944 on 14 and 22 DF,  p-value: 0.002049

lm_features_2 = lm(mean_TRPE ~ xpace + xload + xgender + xage + xheight + xweight + xwaist + xhip + 
                     mean_Wrist_Acc_Std + 
                     mean_Wrist_Jerk_Std +
                     mean_Trunk_Acc_Std + 
                     mean_Trunk_Jerk_Std + 
                     mean_Upperarm_Acc_Std + 
                     mean_Upperarm_Jerk_Std,
    agg_df)

summary(lm_features_2)

Call:
lm(formula = mean_TRPE ~ xpace + xload + xgender + xage + xheight + 
    xweight + xwaist + xhip + mean_Wrist_Acc_Std + mean_Wrist_Jerk_Std + 
    mean_Trunk_Acc_Std + mean_Trunk_Jerk_Std + mean_Upperarm_Acc_Std + 
    mean_Upperarm_Jerk_Std, data = agg_df)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.77857 -0.29064 -0.02397  0.17467  1.29296 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)   
(Intercept)             5.07205    5.85360   0.866  0.39558   
xpace                   0.93933    0.29590   3.175  0.00439 **
xload                  -0.02528    0.05372  -0.471  0.64254   
xgenderF               -2.16410    1.36014  -1.591  0.12586   
xage                    0.09937    0.05850   1.699  0.10348   
xheight                -0.07202    0.03602  -1.999  0.05807 . 
xweight                 0.02038    0.01533   1.329  0.19734   
xwaist                 -0.09331    0.03348  -2.787  0.01074 * 
xhip                    0.09484    0.02527   3.753  0.00110 **
mean_Wrist_Acc_Std     15.44182    6.55016   2.357  0.02771 * 
mean_Wrist_Jerk_Std    -0.36603    0.32257  -1.135  0.26869   
mean_Trunk_Acc_Std     -1.70300   30.18221  -0.056  0.95551   
mean_Trunk_Jerk_Std     0.15461    1.15065   0.134  0.89433   
mean_Upperarm_Acc_Std   6.67274   18.16084   0.367  0.71681   
mean_Upperarm_Jerk_Std -0.24444    0.71293  -0.343  0.73496   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.597 on 22 degrees of freedom
Multiple R-squared:  0.7625,    Adjusted R-squared:  0.6113 
F-statistic: 5.045 on 14 and 22 DF,  p-value: 0.0003881

lm_features_3 = lm(mean_TRPE ~ xpace + xload + xgender + xage + xheight + xweight + xwaist + xhip +
                     mean_Wrist_Acc_Min + 
                     mean_Wrist_Jerk_Min + 
                     mean_Trunk_Acc_Min + 
                     mean_Trunk_Jerk_Min +
                     mean_Upperarm_Acc_Min + 
                     mean_Upperarm_Jerk_Min,
    agg_df)

summary(lm_features_3)

Call:
lm(formula = mean_TRPE ~ xpace + xload + xgender + xage + xheight + 
    xweight + xwaist + xhip + mean_Wrist_Acc_Min + mean_Wrist_Jerk_Min + 
    mean_Trunk_Acc_Min + mean_Trunk_Jerk_Min + mean_Upperarm_Acc_Min + 
    mean_Upperarm_Jerk_Min, data = agg_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.6921 -0.3286 -0.1056  0.3810  1.0867 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)            10.592217   6.368251   1.663  0.11044    
xpace                   0.751068   0.287160   2.616  0.01579 *  
xload                   0.064220   0.053805   1.194  0.24536    
xgenderF               -3.943430   1.312989  -3.003  0.00654 ** 
xage                    0.138528   0.054200   2.556  0.01802 *  
xheight                -0.094804   0.035627  -2.661  0.01427 *  
xweight                 0.006854   0.016401   0.418  0.68007    
xwaist                 -0.126944   0.036394  -3.488  0.00208 ** 
xhip                    0.127879   0.026009   4.917 6.45e-05 ***
mean_Wrist_Acc_Min     -6.482020   2.134654  -3.037  0.00606 ** 
mean_Wrist_Jerk_Min     0.321737  11.776172   0.027  0.97845    
mean_Trunk_Acc_Min      1.921331   2.365321   0.812  0.42533    
mean_Trunk_Jerk_Min    30.308653  27.979542   1.083  0.29043    
mean_Upperarm_Acc_Min   0.206734   1.247112   0.166  0.86985    
mean_Upperarm_Jerk_Min -8.804264  15.792487  -0.557  0.58282    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.5976 on 22 degrees of freedom
Multiple R-squared:  0.7621,    Adjusted R-squared:  0.6107 
F-statistic: 5.033 on 14 and 22 DF,  p-value: 0.0003945

lm_features_4 = lm(mean_TRPE ~ xpace + xload + xgender + xage + xheight + xweight + xwaist + xhip + 
                     mean_Wrist_Acc_Max + 
                     mean_Wrist_Jerk_Max + 
                     mean_Trunk_Acc_Max + 
                     mean_Trunk_Jerk_Max +
                     mean_Upperarm_Acc_Max + 
                     mean_Upperarm_Jerk_Max,
    agg_df)

summary(lm_features_4)

Call:
lm(formula = mean_TRPE ~ xpace + xload + xgender + xage + xheight + 
    xweight + xwaist + xhip + mean_Wrist_Acc_Max + mean_Wrist_Jerk_Max + 
    mean_Trunk_Acc_Max + mean_Trunk_Jerk_Max + mean_Upperarm_Acc_Max + 
    mean_Upperarm_Jerk_Max, data = agg_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.6753 -0.3152 -0.1662  0.3223  1.2635 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)   
(Intercept)             4.53037    7.35331   0.616  0.54415   
xpace                   0.85971    0.31540   2.726  0.01234 * 
xload                   0.09258    0.03232   2.865  0.00901 **
xgenderF               -2.17886    1.16024  -1.878  0.07371 . 
xage                    0.11017    0.06247   1.764  0.09167 . 
xheight                -0.07195    0.03656  -1.968  0.06179 . 
xweight                 0.03443    0.01756   1.961  0.06270 . 
xwaist                 -0.10685    0.03533  -3.024  0.00623 **
xhip                    0.08173    0.02548   3.208  0.00406 **
mean_Wrist_Acc_Max      0.89891    0.71827   1.251  0.22390   
mean_Wrist_Jerk_Max    -0.01132    0.02907  -0.389  0.70077   
mean_Trunk_Acc_Max     -0.66051    2.83164  -0.233  0.81772   
mean_Trunk_Jerk_Max    -0.01779    0.08611  -0.207  0.83824   
mean_Upperarm_Acc_Max   2.38047    2.29483   1.037  0.31086   
mean_Upperarm_Jerk_Max -0.08223    0.06181  -1.330  0.19702   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6272 on 22 degrees of freedom
Multiple R-squared:  0.7379,    Adjusted R-squared:  0.571 
F-statistic: 4.423 on 14 and 22 DF,  p-value: 0.0009662

lm_features_5 = lm(mean_TRPE ~ xpace + xload + xgender + xage + xheight + xweight + xwaist + xhip + 
                     mean_Wrist_Acc_Med + 
                     mean_Wrist_Jerk_Med + 
                     mean_Trunk_Acc_Med + 
                     mean_Trunk_Jerk_Med +
                     mean_Upperarm_Acc_Med + 
                     mean_Upperarm_Jerk_Med,
    agg_df)

summary(lm_features_5)

Call:
lm(formula = mean_TRPE ~ xpace + xload + xgender + xage + xheight + 
    xweight + xwaist + xhip + mean_Wrist_Acc_Med + mean_Wrist_Jerk_Med + 
    mean_Trunk_Acc_Med + mean_Trunk_Jerk_Med + mean_Upperarm_Acc_Med + 
    mean_Upperarm_Jerk_Med, data = agg_df)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.70202 -0.28089 -0.09378  0.36267  1.06421 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)            -52.34293   27.24784  -1.921 0.067785 .  
xpace                    0.93710    0.29466   3.180 0.004329 ** 
xload                   -0.04543    0.07275  -0.625 0.538708    
xgenderF                -1.50918    1.01432  -1.488 0.150979    
xage                     0.07891    0.05308   1.487 0.151326    
xheight                 -0.05066    0.02851  -1.777 0.089407 .  
xweight                  0.01275    0.01603   0.795 0.435103    
xwaist                  -0.05317    0.02910  -1.827 0.081234 .  
xhip                     0.09389    0.02245   4.183 0.000386 ***
mean_Wrist_Acc_Med      48.37228   20.17770   2.397 0.025443 *  
mean_Wrist_Jerk_Med      0.18196    0.34376   0.529 0.601883    
mean_Trunk_Acc_Med      15.52905   10.33003   1.503 0.146984    
mean_Trunk_Jerk_Med      0.54350    0.64604   0.841 0.409244    
mean_Upperarm_Acc_Med  -11.33600   10.70835  -1.059 0.301261    
mean_Upperarm_Jerk_Med  -0.03150    0.43812  -0.072 0.943334    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6016 on 22 degrees of freedom
Multiple R-squared:  0.7588,    Adjusted R-squared:  0.6053 
F-statistic: 4.944 on 14 and 22 DF,  p-value: 0.000448

lm_features_6 = lm(mean_TRPE ~ xpace + xload + xgender + xage + xheight + xweight + xwaist + xhip + 
                     mean_Wrist_Acc_Mad + 
                     mean_Wrist_Jerk_Mad + 
                     mean_Trunk_Acc_Mad + 
                     mean_Trunk_Jerk_Mad +
                     mean_Upperarm_Acc_Mad + 
                     mean_Upperarm_Jerk_Mad,
    agg_df)

summary(lm_features_6)

Call:
lm(formula = mean_TRPE ~ xpace + xload + xgender + xage + xheight + 
    xweight + xwaist + xhip + mean_Wrist_Acc_Mad + mean_Wrist_Jerk_Mad + 
    mean_Trunk_Acc_Mad + mean_Trunk_Jerk_Mad + mean_Upperarm_Acc_Mad + 
    mean_Upperarm_Jerk_Mad, data = agg_df)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.67918 -0.35736 -0.00721  0.17209  1.51644 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)   
(Intercept)              3.45436    5.62339   0.614  0.54533   
xpace                    0.86734    0.29399   2.950  0.00740 **
xload                   -0.08200    0.07840  -1.046  0.30696   
xgenderF                -1.99626    1.30901  -1.525  0.14151   
xage                     0.09713    0.05929   1.638  0.11561   
xheight                 -0.06445    0.03518  -1.832  0.08049 . 
xweight                  0.01749    0.01561   1.120  0.27467   
xwaist                  -0.07706    0.03346  -2.303  0.03112 * 
xhip                     0.09103    0.02609   3.490  0.00207 **
mean_Wrist_Acc_Mad      29.89322   12.36910   2.417  0.02440 * 
mean_Wrist_Jerk_Mad     -0.77317    0.56356  -1.372  0.18392   
mean_Trunk_Acc_Mad       5.38496   33.99600   0.158  0.87559   
mean_Trunk_Jerk_Mad      0.16636    1.23470   0.135  0.89404   
mean_Upperarm_Acc_Mad  -25.91824   30.91308  -0.838  0.41081   
mean_Upperarm_Jerk_Mad   1.10837    1.18773   0.933  0.36086   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6145 on 22 degrees of freedom
Multiple R-squared:  0.7484,    Adjusted R-squared:  0.5883 
F-statistic: 4.674 on 14 and 22 DF,  p-value: 0.0006635

lm_features_7 = lm(mean_TRPE ~ xpace + xload + xgender + xage + xheight + xweight + xwaist + xhip + 
                     mean_Wrist_Acc_Rms + 
                     mean_Wrist_Jerk_Rms + 
                     mean_Trunk_Acc_Rms + 
                     mean_Trunk_Jerk_Rms +
                     mean_Upperarm_Acc_Rms + 
                     mean_Upperarm_Jerk_Rms,
    agg_df)

summary(lm_features_7)

Call:
lm(formula = mean_TRPE ~ xpace + xload + xgender + xage + xheight + 
    xweight + xwaist + xhip + mean_Wrist_Acc_Rms + mean_Wrist_Jerk_Rms + 
    mean_Trunk_Acc_Rms + mean_Trunk_Jerk_Rms + mean_Upperarm_Acc_Rms + 
    mean_Upperarm_Jerk_Rms, data = agg_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.6710 -0.3429 -0.1620  0.3176  1.0628 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)            -36.32449   16.72335  -2.172 0.040904 *  
xpace                    1.09655    0.28613   3.832 0.000907 ***
xload                    0.05790    0.03986   1.453 0.160445    
xgenderF                -3.56263    1.09625  -3.250 0.003673 ** 
xage                     0.14808    0.05296   2.796 0.010528 *  
xheight                 -0.11620    0.03258  -3.567 0.001725 ** 
xweight                  0.02167    0.01688   1.284 0.212382    
xwaist                  -0.12171    0.03084  -3.946 0.000687 ***
xhip                     0.12149    0.02346   5.178 3.43e-05 ***
mean_Wrist_Acc_Rms      41.27947   13.70478   3.012 0.006413 ** 
mean_Wrist_Jerk_Rms     -0.45704    0.28084  -1.627 0.117885    
mean_Trunk_Acc_Rms      16.06923    9.60972   1.672 0.108652    
mean_Trunk_Jerk_Rms      0.05327    0.31416   0.170 0.866916    
mean_Upperarm_Acc_Rms   -7.96835    9.45170  -0.843 0.408268    
mean_Upperarm_Jerk_Rms  -0.02579    0.23472  -0.110 0.913518    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.5754 on 22 degrees of freedom
Multiple R-squared:  0.7794,    Adjusted R-squared:  0.639 
F-statistic: 5.551 on 14 and 22 DF,  p-value: 0.0001938

lm_features_8 = lm(mean_TRPE ~ xpace + xload + xgender + xage + xheight + xweight + xwaist + xhip + 
                     mean_Wrist_Acc_Perc25 + 
                     mean_Wrist_Jerk_Perc25 + 
                     mean_Trunk_Acc_Perc25 + 
                     mean_Trunk_Jerk_Perc25 +
                     mean_Upperarm_Acc_Perc25 + 
                     mean_Upperarm_Jerk_Perc25,
    agg_df)

summary(lm_features_8)

Call:
lm(formula = mean_TRPE ~ xpace + xload + xgender + xage + xheight + 
    xweight + xwaist + xhip + mean_Wrist_Acc_Perc25 + mean_Wrist_Jerk_Perc25 + 
    mean_Trunk_Acc_Perc25 + mean_Trunk_Jerk_Perc25 + mean_Upperarm_Acc_Perc25 + 
    mean_Upperarm_Jerk_Perc25, data = agg_df)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.84449 -0.36499 -0.07673  0.31656  1.56710 

Coefficients:
                            Estimate Std. Error t value Pr(>|t|)    
(Intercept)                -2.023162  17.086175  -0.118 0.906818    
xpace                       0.793563   0.308290   2.574 0.017310 *  
xload                       0.010602   0.086765   0.122 0.903860    
xgenderF                   -2.137988   1.254368  -1.704 0.102384    
xage                        0.110812   0.058013   1.910 0.069238 .  
xheight                    -0.046546   0.034385  -1.354 0.189586    
xweight                     0.009553   0.019523   0.489 0.629456    
xwaist                     -0.068556   0.033199  -2.065 0.050907 .  
xhip                        0.103005   0.026947   3.823 0.000929 ***
mean_Wrist_Acc_Perc25       8.341082   8.810523   0.947 0.354066    
mean_Wrist_Jerk_Perc25      0.491991   0.674914   0.729 0.473710    
mean_Trunk_Acc_Perc25      10.125576  12.090900   0.837 0.411343    
mean_Trunk_Jerk_Perc25      2.537041   1.689372   1.502 0.147376    
mean_Upperarm_Acc_Perc25  -17.226105  10.880735  -1.583 0.127652    
mean_Upperarm_Jerk_Perc25  -0.212549   1.109622  -0.192 0.849851    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6452 on 22 degrees of freedom
Multiple R-squared:  0.7226,    Adjusted R-squared:  0.5461 
F-statistic: 4.094 on 14 and 22 DF,  p-value: 0.001612

lm_features_9 = lm(mean_TRPE ~ xpace + xload + xgender + xage + xheight + xweight + xwaist + xhip + 
                     mean_Wrist_Acc_Perc75 + 
                     mean_Wrist_Jerk_Perc75 + 
                     mean_Trunk_Acc_Perc75 + 
                     mean_Trunk_Jerk_Perc75 +
                     mean_Upperarm_Acc_Perc75 + 
                     mean_Upperarm_Jerk_Perc75,
    agg_df)

summary(lm_features_9)

Call:
lm(formula = mean_TRPE ~ xpace + xload + xgender + xage + xheight + 
    xweight + xwaist + xhip + mean_Wrist_Acc_Perc75 + mean_Wrist_Jerk_Perc75 + 
    mean_Trunk_Acc_Perc75 + mean_Trunk_Jerk_Perc75 + mean_Upperarm_Acc_Perc75 + 
    mean_Upperarm_Jerk_Perc75, data = agg_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.8573 -0.2734 -0.1084  0.2556  1.2298 

Coefficients:
                            Estimate Std. Error t value Pr(>|t|)   
(Intercept)               -21.791649  16.154048  -1.349  0.19106   
xpace                       0.971966   0.298717   3.254  0.00364 **
xload                      -0.096535   0.088574  -1.090  0.28756   
xgenderF                   -1.197924   1.160977  -1.032  0.31336   
xage                        0.100537   0.055650   1.807  0.08452 . 
xheight                    -0.041237   0.031065  -1.327  0.19798   
xweight                     0.005523   0.016423   0.336  0.73986   
xwaist                     -0.039343   0.032852  -1.198  0.24382   
xhip                        0.082618   0.024780   3.334  0.00301 **
mean_Wrist_Acc_Perc75      23.674933   9.044782   2.618  0.01572 * 
mean_Wrist_Jerk_Perc75     -0.331407   0.263287  -1.259  0.22132   
mean_Trunk_Acc_Perc75       9.895349   9.213907   1.074  0.29448   
mean_Trunk_Jerk_Perc75      0.367945   0.374583   0.982  0.33664   
mean_Upperarm_Acc_Perc75  -13.476362  10.335148  -1.304  0.20574   
mean_Upperarm_Jerk_Perc75   0.397133   0.330956   1.200  0.24292   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6084 on 22 degrees of freedom
Multiple R-squared:  0.7533,    Adjusted R-squared:  0.5964 
F-statistic: 4.799 on 14 and 22 DF,  p-value: 0.0005518

lm_features_10 = lm(mean_TRPE ~ xpace + xload + xgender + xage + xheight + xweight + xwaist + xhip +
                     mean_Wrist_Acc_Min + 
                     mean_Wrist_Jerk_Min + 
                     mean_Trunk_Acc_Min + 
                     mean_Trunk_Jerk_Min +
                     mean_Upperarm_Acc_Min + 
                     mean_Upperarm_Jerk_Min +
                     mean_Wrist_Acc_Max + 
                     mean_Wrist_Jerk_Max + 
                     mean_Trunk_Acc_Max + 
                     mean_Trunk_Jerk_Max +
                     mean_Upperarm_Acc_Max + 
                     mean_Upperarm_Jerk_Max,
                      
    agg_df)

summary(lm_features_10)

Call:
lm(formula = mean_TRPE ~ xpace + xload + xgender + xage + xheight + 
    xweight + xwaist + xhip + mean_Wrist_Acc_Min + mean_Wrist_Jerk_Min + 
    mean_Trunk_Acc_Min + mean_Trunk_Jerk_Min + mean_Upperarm_Acc_Min + 
    mean_Upperarm_Jerk_Min + mean_Wrist_Acc_Max + mean_Wrist_Jerk_Max + 
    mean_Trunk_Acc_Max + mean_Trunk_Jerk_Max + mean_Upperarm_Acc_Max + 
    mean_Upperarm_Jerk_Max, data = agg_df)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.60570 -0.35291 -0.09057  0.29447  1.06320 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)   
(Intercept)             7.6962351 12.9701433   0.593  0.56122   
xpace                   0.8211095  0.3709158   2.214  0.04172 * 
xload                   0.0438230  0.0748629   0.585  0.56646   
xgenderF               -3.9066886  1.9104134  -2.045  0.05767 . 
xage                    0.1407816  0.0930590   1.513  0.14983   
xheight                -0.0944751  0.0544851  -1.734  0.10215   
xweight                 0.0195332  0.0288560   0.677  0.50813   
xwaist                 -0.1387497  0.0492230  -2.819  0.01235 * 
xhip                    0.1228654  0.0391383   3.139  0.00634 **
mean_Wrist_Acc_Min     -4.5143489  4.3758647  -1.032  0.31758   
mean_Wrist_Jerk_Min    -1.9981976 14.8806342  -0.134  0.89485   
mean_Trunk_Acc_Min      1.6733716  7.0638807   0.237  0.81575   
mean_Trunk_Jerk_Min    46.8198676 41.1729108   1.137  0.27222   
mean_Upperarm_Acc_Min   1.1179984  2.1872858   0.511  0.61624   
mean_Upperarm_Jerk_Min -8.8318869 21.7987187  -0.405  0.69073   
mean_Wrist_Acc_Max      0.2984092  1.3093894   0.228  0.82261   
mean_Wrist_Jerk_Max    -0.0001444  0.0390171  -0.004  0.99709   
mean_Trunk_Acc_Max      0.3557997  3.4512628   0.103  0.91917   
mean_Trunk_Jerk_Max    -0.0332983  0.1163733  -0.286  0.77845   
mean_Upperarm_Acc_Max   1.1194518  2.9809839   0.376  0.71220   
mean_Upperarm_Jerk_Max -0.0153030  0.0962198  -0.159  0.87563   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6776 on 16 degrees of freedom
Multiple R-squared:  0.7775,    Adjusted R-squared:  0.4994 
F-statistic: 2.795 on 20 and 16 DF,  p-value: 0.02071

lm_features_11 = lm(mean_TRPE ~ xpace + xload + xgender + xage + xheight + xweight + xwaist + xhip + 
                     mean_Wrist_Acc_Perc25 + 
                     mean_Wrist_Jerk_Perc25 + 
                     mean_Trunk_Acc_Perc25 + 
                     mean_Trunk_Jerk_Perc25 +
                     mean_Upperarm_Acc_Perc25 + 
                     mean_Upperarm_Jerk_Perc25 +
                     mean_Wrist_Acc_Perc75 + 
                     mean_Wrist_Jerk_Perc75 + 
                     mean_Trunk_Acc_Perc75 + 
                     mean_Trunk_Jerk_Perc75 +
                     mean_Upperarm_Acc_Perc75 + 
                     mean_Upperarm_Jerk_Perc75,
    agg_df)

summary(lm_features_11)

Call:
lm(formula = mean_TRPE ~ xpace + xload + xgender + xage + xheight + 
    xweight + xwaist + xhip + mean_Wrist_Acc_Perc25 + mean_Wrist_Jerk_Perc25 + 
    mean_Trunk_Acc_Perc25 + mean_Trunk_Jerk_Perc25 + mean_Upperarm_Acc_Perc25 + 
    mean_Upperarm_Jerk_Perc25 + mean_Wrist_Acc_Perc75 + mean_Wrist_Jerk_Perc75 + 
    mean_Trunk_Acc_Perc75 + mean_Trunk_Jerk_Perc75 + mean_Upperarm_Acc_Perc75 + 
    mean_Upperarm_Jerk_Perc75, data = agg_df)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.85811 -0.23054 -0.03207  0.17146  1.23273 

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)   
(Intercept)               -36.44448   25.27124  -1.442  0.16855   
xpace                       1.12522    0.33260   3.383  0.00379 **
xload                      -0.14472    0.12341  -1.173  0.25811   
xgenderF                   -1.68809    1.37044  -1.232  0.23582   
xage                        0.13852    0.06477   2.139  0.04825 * 
xheight                    -0.04943    0.03680  -1.343  0.19800   
xweight                     0.01212    0.02033   0.596  0.55923   
xwaist                     -0.04155    0.03613  -1.150  0.26701   
xhip                        0.08011    0.02880   2.781  0.01335 * 
mean_Wrist_Acc_Perc25       5.82471   12.41271   0.469  0.64522   
mean_Wrist_Jerk_Perc25     -0.17871    1.19639  -0.149  0.88312   
mean_Trunk_Acc_Perc25      84.57578   51.03672   1.657  0.11696   
mean_Trunk_Jerk_Perc25      9.49745    7.41082   1.282  0.21826   
mean_Upperarm_Acc_Perc25    3.48359   29.56851   0.118  0.90768   
mean_Upperarm_Jerk_Perc25  -0.84426    2.61219  -0.323  0.75073   
mean_Wrist_Acc_Perc75      32.94292   12.51798   2.632  0.01814 * 
mean_Wrist_Jerk_Perc75     -0.41969    0.32527  -1.290  0.21529   
mean_Trunk_Acc_Perc75     -65.85007   48.84353  -1.348  0.19638   
mean_Trunk_Jerk_Perc75      0.26739    0.69393   0.385  0.70507   
mean_Upperarm_Acc_Perc75  -26.22405   34.49721  -0.760  0.45820   
mean_Upperarm_Jerk_Perc75   0.55594    0.56031   0.992  0.33587   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6258 on 16 degrees of freedom
Multiple R-squared:  0.8102,    Adjusted R-squared:  0.5729 
F-statistic: 3.415 on 20 and 16 DF,  p-value: 0.007918

5 Plotting transformed RPE and sensors’ median

6 Plotting the decimated signals

7 Plotting a slice of decimated signals for subject 9 and condition 2.5-15

8 Check the correlation between the scalar inputs and the PCs from the 30-second intervals

# The canonical correlation using cca.fd::fda

smooth_functional_var = function(nbasis, var){
  ncol=dim(var)[2]
  basis_var <- fda::create.bspline.basis(c(0,ncol), nbasis= nbasis)
  smooth_var <- fda::smooth.basis(argvals = seq(0,ncol, length.out = dim(t(var))[1]), y = t(var), fdParobj = basis_var)$fd
  list = list(basis = basis_var,
              smooth_var = smooth_var)
  return(list)
  }

# change the format of acceleration PCs to a fd object
trunk_AccX_fd = smooth_functional_var(dim(input_data_fpca_30sec$trunk_Accelerometer_X)[2],input_data_fpca_30sec$trunk_Accelerometer_X)$smooth_var
trunk_AccY_fd = smooth_functional_var(dim(input_data_fpca_30sec$trunk_Accelerometer_Y)[2],input_data_fpca_30sec$trunk_Accelerometer_Y)$smooth_var
trunk_AccZ_fd = smooth_functional_var(dim(input_data_fpca_30sec$trunk_Accelerometer_Z)[2],input_data_fpca_30sec$trunk_Accelerometer_Z)$smooth_var
wrist_AccX_fd = smooth_functional_var(dim(input_data_fpca_30sec$wrist_Accelerometer_X)[2],input_data_fpca_30sec$wrist_Accelerometer_X)$smooth_var
wrist_AccY_fd = smooth_functional_var(dim(input_data_fpca_30sec$wrist_Accelerometer_Y)[2],input_data_fpca_30sec$wrist_Accelerometer_Y)$smooth_var
wrist_AccZ_fd = smooth_functional_var(dim(input_data_fpca_30sec$wrist_Accelerometer_Z)[2],input_data_fpca_30sec$wrist_Accelerometer_Z)$smooth_var
upperarm_AccX_fd = smooth_functional_var(dim(input_data_fpca_30sec$upperarm_Accelerometer_X)[2],input_data_fpca_30sec$upperarm_Accelerometer_X)$smooth_var
upperarm_AccY_fd = smooth_functional_var(dim(input_data_fpca_30sec$upperarm_Accelerometer_Y)[2],input_data_fpca_30sec$upperarm_Accelerometer_Y)$smooth_var
upperarm_AccZ_fd = smooth_functional_var(dim(input_data_fpca_30sec$upperarm_Accelerometer_Z)[2],input_data_fpca_30sec$upperarm_Accelerometer_Z)$smooth_var

# change the format of scalar inputs to an fd object (considering replicates of them over 4 points)

load_fd = smooth_functional_var(4,matrix(data=rep(input_data_fpca_30sec$xload,4),nrow=44,ncol=4))$smooth_var
pace_fd = smooth_functional_var(4,matrix(data=rep(input_data_fpca_30sec$xpace,4),nrow=44,ncol=4))$smooth_var
height_fd = smooth_functional_var(4,matrix(data=rep(input_data_fpca_30sec$xheight,4),nrow=44,ncol=4))$smooth_var
weight_fd = smooth_functional_var(4,matrix(data=rep(input_data_fpca_30sec$xweight,4),nrow=44,ncol=4))$smooth_var
age_fd = smooth_functional_var(4,matrix(data=rep(input_data_fpca_30sec$xage,4),nrow=44,ncol=4))$smooth_var
hip_fd = smooth_functional_var(4,matrix(data=rep(input_data_fpca_30sec$xhip,4),nrow=44,ncol=4))$smooth_var
waist_fd = smooth_functional_var(4,matrix(data=rep(input_data_fpca_30sec$xwaist,4),nrow=44,ncol=4))$smooth_var
# gender couldn't be changed to a functional variable as it is not numeric. May be we can replace the factors with numbers?

# finding the correlation using the fda::cor.fd

##for Trunk_acc_x

  
  corfd_load_trunkACCX <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$trunk_Accelerometer_X)[2], length=dim(input_data_fpca_30sec$trunk_Accelerometer_X)[2]), trunk_AccX_fd,
  seq(0, 4, length=4), load_fd)[,1]
  corfd_pace_trunkACCX <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$trunk_Accelerometer_X)[2], length=dim(input_data_fpca_30sec$trunk_Accelerometer_X)[2]), trunk_AccX_fd,
  seq(0, 4, length=4), pace_fd)[,1]
  corfd_height_trunkACCX <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$trunk_Accelerometer_X)[2], length=dim(input_data_fpca_30sec$trunk_Accelerometer_X)[2]), trunk_AccX_fd,
  seq(0, 4, length=4), height_fd)[,1]
  corfd_weight_trunkACCX <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$trunk_Accelerometer_X)[2], length=dim(input_data_fpca_30sec$trunk_Accelerometer_X)[2]), trunk_AccX_fd,
  seq(0, 4, length=4), weight_fd)[,1]
  corfd_hip_trunkACCX <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$trunk_Accelerometer_X)[2], length=dim(input_data_fpca_30sec$trunk_Accelerometer_X)[2]), trunk_AccX_fd,
  seq(0, 4, length=4), hip_fd)[,1]
  corfd_waist_trunkACCX <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$trunk_Accelerometer_X)[2], length=dim(input_data_fpca_30sec$trunk_Accelerometer_X)[2]), trunk_AccX_fd,
  seq(0, 4, length=4), waist_fd)[,1]
  PC_no = seq(1,dim(input_data_fpca_30sec$trunk_Accelerometer_X)[2],1)
  trunkACCX_corfd_summary = data.frame(PC_no,corfd_load_trunkACCX, corfd_pace_trunkACCX, corfd_height_trunkACCX, corfd_weight_trunkACCX, corfd_hip_trunkACCX, corfd_waist_trunkACCX )
  trunkACCX_corfd_summary_long = pivot_longer(trunkACCX_corfd_summary, cols = contains('trunkACCX') ,names_to = 'factors', values_to = 'corr_value')
  
  
  trunkACCX_corfd_plot = ggplot(data=trunkACCX_corfd_summary_long, aes(x= PC_no, y= corr_value))+
  geom_line() +
  geom_point()+ 
    facet_wrap(~factors, ncol =3)+
  ggtitle('Correlation coefficients of different factors and PCs of trunk ACC on X direction')+ 
  ggplot2::labs(
    x = 'PCs', y = 'Correlation value') +
  ggplot2::theme(legend.position = 'bottom') +
  ggplot2::theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  )+
  ggplot2::theme(
     axis.title = ggplot2::element_text(size = 14),  # Increase axis label size
     axis.text = ggplot2::element_text(size = 14) # Increase axis text size
  )
  
 print(trunkACCX_corfd_plot)

  ##for Trunk_acc_y
  
  corfd_load_trunkACCY <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$trunk_Accelerometer_Y)[2], length=dim(input_data_fpca_30sec$trunk_Accelerometer_Y)[2]), trunk_AccY_fd,
  seq(0, 4, length=4), load_fd)[,1]
  corfd_pace_trunkACCY <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$trunk_Accelerometer_Y)[2], length=dim(input_data_fpca_30sec$trunk_Accelerometer_Y)[2]), trunk_AccY_fd,
  seq(0, 4, length=4), pace_fd)[,1]
  corfd_height_trunkACCY <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$trunk_Accelerometer_Y)[2], length=dim(input_data_fpca_30sec$trunk_Accelerometer_Y)[2]), trunk_AccY_fd,
  seq(0, 4, length=4), height_fd)[,1]
  corfd_weight_trunkACCY <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$trunk_Accelerometer_Y)[2], length=dim(input_data_fpca_30sec$trunk_Accelerometer_Y)[2]), trunk_AccY_fd,
  seq(0, 4, length=4), weight_fd)[,1]
  corfd_hip_trunkACCY <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$trunk_Accelerometer_Y)[2], length=dim(input_data_fpca_30sec$trunk_Accelerometer_Y)[2]), trunk_AccY_fd,
  seq(0, 4, length=4), hip_fd)[,1]
  corfd_waist_trunkACCY <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$trunk_Accelerometer_Y)[2], length=dim(input_data_fpca_30sec$trunk_Accelerometer_Y)[2]), trunk_AccY_fd,
  seq(0, 4, length=4), waist_fd)[,1]
  PC_no = seq(1,dim(input_data_fpca_30sec$trunk_Accelerometer_Y)[2],1)
  trunkACCY_corfd_summary = data.frame(PC_no,corfd_load_trunkACCY, corfd_pace_trunkACCY, corfd_height_trunkACCY, corfd_weight_trunkACCY, corfd_hip_trunkACCY, corfd_waist_trunkACCY )
  trunkACCY_corfd_summary_long = pivot_longer(trunkACCY_corfd_summary, cols = contains('trunkACCY') ,names_to = 'factors', values_to = 'corr_value')
  
  
  trunkACCY_corfd_plot = ggplot(data=trunkACCY_corfd_summary_long, aes(x= PC_no, y= corr_value))+
  geom_line() +
  geom_point()+ 
    facet_wrap(~factors, ncol =3)+
  ggtitle('Correlation coefficients of different factors and PCs of trunk ACC on Y direction')+ 
  ggplot2::labs(
    x = 'PCs', y = 'Correlation value') +
  ggplot2::theme(legend.position = 'bottom') +
  ggplot2::theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  )+
  ggplot2::theme(
     axis.title = ggplot2::element_text(size = 14),  # Increase axis label size
     axis.text = ggplot2::element_text(size = 14) # Increase axis text size
  )
  
  print(trunkACCY_corfd_plot)

 ##for Trunk_acc_z
  
  corfd_load_trunkACCZ <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$trunk_Accelerometer_Z)[2], length=dim(input_data_fpca_30sec$trunk_Accelerometer_Z)[2]), trunk_AccZ_fd,
  seq(0, 4, length=4), load_fd)[,1]
  corfd_pace_trunkACCZ <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$trunk_Accelerometer_Z)[2], length=dim(input_data_fpca_30sec$trunk_Accelerometer_Z)[2]), trunk_AccZ_fd,
  seq(0, 4, length=4), pace_fd)[,1]
  corfd_height_trunkACCZ <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$trunk_Accelerometer_Z)[2], length=dim(input_data_fpca_30sec$trunk_Accelerometer_Z)[2]), trunk_AccZ_fd,
  seq(0, 4, length=4), height_fd)[,1]
  corfd_weight_trunkACCZ <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$trunk_Accelerometer_Z)[2], length=dim(input_data_fpca_30sec$trunk_Accelerometer_Z)[2]), trunk_AccZ_fd,
  seq(0, 4, length=4), weight_fd)[,1]
  corfd_hip_trunkACCZ <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$trunk_Accelerometer_Z)[2], length=dim(input_data_fpca_30sec$trunk_Accelerometer_Z)[2]), trunk_AccZ_fd,
  seq(0, 4, length=4), hip_fd)[,1]
  corfd_waist_trunkACCZ <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$trunk_Accelerometer_Z)[2], length=dim(input_data_fpca_30sec$trunk_Accelerometer_Z)[2]), trunk_AccZ_fd,
  seq(0, 4, length=4), waist_fd)[,1]
  PC_no = seq(1,dim(input_data_fpca_30sec$trunk_Accelerometer_Z)[2],1)
  trunkACCZ_corfd_summary = data.frame(PC_no,corfd_load_trunkACCZ, corfd_pace_trunkACCZ, corfd_height_trunkACCZ, corfd_weight_trunkACCZ, corfd_hip_trunkACCZ, corfd_waist_trunkACCZ )
  trunkACCZ_corfd_summary_long = pivot_longer(trunkACCZ_corfd_summary, cols = contains('trunkACCZ') ,names_to = 'factors', values_to = 'corr_value')
  
  
  trunkACCZ_corfd_plot = ggplot(data=trunkACCZ_corfd_summary_long, aes(x= PC_no, y= corr_value))+
  geom_line() +
  geom_point()+ 
    facet_wrap(~factors, ncol =3)+
  ggtitle('Correlation coefficients of different factors and PCs of trunk ACC on Z direction')+ 
  ggplot2::labs(
    x = 'PCs', y = 'Correlation value') +
  ggplot2::theme(legend.position = 'bottom') +
  ggplot2::theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  )+
  ggplot2::theme(
     axis.title = ggplot2::element_text(size = 14),  # Increase axis label size
     axis.text = ggplot2::element_text(size = 14) # Increase axis text size
  )
  
  print(trunkACCZ_corfd_plot)

  ##for Wrist_acc_x
  
  corfd_load_wristACCX <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$wrist_Accelerometer_X)[2], length=dim(input_data_fpca_30sec$wrist_Accelerometer_X)[2]), wrist_AccX_fd,
  seq(0, 4, length=4), load_fd)[,1]
  corfd_pace_wristACCX <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$wrist_Accelerometer_X)[2], length=dim(input_data_fpca_30sec$wrist_Accelerometer_X)[2]), wrist_AccX_fd,
  seq(0, 4, length=4), pace_fd)[,1]
  corfd_height_wristACCX <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$wrist_Accelerometer_X)[2], length=dim(input_data_fpca_30sec$wrist_Accelerometer_X)[2]), wrist_AccX_fd,
  seq(0, 4, length=4), height_fd)[,1]
  corfd_weight_wristACCX <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$wrist_Accelerometer_X)[2], length=dim(input_data_fpca_30sec$wrist_Accelerometer_X)[2]), wrist_AccX_fd,
  seq(0, 4, length=4), weight_fd)[,1]
  corfd_hip_wristACCX <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$wrist_Accelerometer_X)[2], length=dim(input_data_fpca_30sec$wrist_Accelerometer_X)[2]), wrist_AccX_fd,
  seq(0, 4, length=4), hip_fd)[,1]
  corfd_waist_wristACCX <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$wrist_Accelerometer_X)[2], length=dim(input_data_fpca_30sec$wrist_Accelerometer_X)[2]), wrist_AccX_fd,
  seq(0, 4, length=4), waist_fd)[,1]
  PC_no = seq(1,dim(input_data_fpca_30sec$wrist_Accelerometer_X)[2],1)
  wristACCX_corfd_summary = data.frame(PC_no,corfd_load_wristACCX, corfd_pace_wristACCX, corfd_height_wristACCX, corfd_weight_wristACCX, corfd_hip_wristACCX, corfd_waist_wristACCX )
  wristACCX_corfd_summary_long = pivot_longer(wristACCX_corfd_summary, cols = contains('wristACCX') ,names_to = 'factors', values_to = 'corr_value')
  
  wristACCX_corfd_plot = ggplot(data=wristACCX_corfd_summary_long, aes(x= PC_no, y= corr_value))+
  geom_line() +
  geom_point()+ 
    facet_wrap(~factors, ncol =3)+
  ggtitle('Correlation coefficients of different factors and PCs of wrist ACC on X direction')+ 
  ggplot2::labs(
    x = 'PCs', y = 'Correlation value') +
  ggplot2::theme(legend.position = 'bottom') +
  ggplot2::theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  )+
  ggplot2::theme(
     axis.title = ggplot2::element_text(size = 14),  # Increase axis label size
     axis.text = ggplot2::element_text(size = 14) # Increase axis text size
  )
  print(wristACCX_corfd_plot)

  ##for Wrist_acc_y
  
  corfd_load_wristACCY <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$wrist_Accelerometer_Y)[2], length=dim(input_data_fpca_30sec$wrist_Accelerometer_Y)[2]), wrist_AccY_fd,
  seq(0, 4, length=4), load_fd)[,1]
  corfd_pace_wristACCY <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$wrist_Accelerometer_Y)[2], length=dim(input_data_fpca_30sec$wrist_Accelerometer_Y)[2]), wrist_AccY_fd,
  seq(0, 4, length=4), pace_fd)[,1]
  corfd_height_wristACCY <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$wrist_Accelerometer_Y)[2], length=dim(input_data_fpca_30sec$wrist_Accelerometer_Y)[2]), wrist_AccY_fd,
  seq(0, 4, length=4), height_fd)[,1]
  corfd_weight_wristACCY <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$wrist_Accelerometer_Y)[2], length=dim(input_data_fpca_30sec$wrist_Accelerometer_Y)[2]), wrist_AccY_fd,
  seq(0, 4, length=4), weight_fd)[,1]
  corfd_hip_wristACCY <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$wrist_Accelerometer_Y)[2], length=dim(input_data_fpca_30sec$wrist_Accelerometer_Y)[2]), wrist_AccY_fd,
  seq(0, 4, length=4), hip_fd)[,1]
  corfd_waist_wristACCY <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$wrist_Accelerometer_Y)[2], length=dim(input_data_fpca_30sec$wrist_Accelerometer_Y)[2]), wrist_AccY_fd,
  seq(0, 4, length=4), waist_fd)[,1]
  PC_no = seq(1,dim(input_data_fpca_30sec$wrist_Accelerometer_Y)[2],1)
  wristACCY_corfd_summary = data.frame(PC_no,corfd_load_wristACCY, corfd_pace_wristACCY, corfd_height_wristACCY, corfd_weight_wristACCY, corfd_hip_wristACCY, corfd_waist_wristACCY )
  wristACCY_corfd_summary_long = pivot_longer(wristACCY_corfd_summary, cols = contains('wristACCY') ,names_to = 'factors', values_to = 'corr_value')
  
  wristACCY_corfd_plot = ggplot(data=wristACCY_corfd_summary_long, aes(x= PC_no, y= corr_value))+
  geom_line() +
  geom_point()+ 
    facet_wrap(~factors, ncol =3)+
  ggtitle('Correlation coefficients of different factors and PCs of wrist ACC on Y direction')+ 
  ggplot2::labs(
    x = 'PCs', y = 'Correlation value') +
  ggplot2::theme(legend.position = 'bottom') +
  ggplot2::theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  )+
  ggplot2::theme(
     axis.title = ggplot2::element_text(size = 14),  # Increase axis label size
     axis.text = ggplot2::element_text(size = 14) # Increase axis text size
  )
  
   print(wristACCY_corfd_plot)

 ##for Wrist_acc_z
  
  corfd_load_wristACCZ <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$wrist_Accelerometer_Z)[2], length=dim(input_data_fpca_30sec$wrist_Accelerometer_Z)[2]), wrist_AccZ_fd,
  seq(0, 4, length=4), load_fd)[,1]
  corfd_pace_wristACCZ <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$wrist_Accelerometer_Z)[2], length=dim(input_data_fpca_30sec$wrist_Accelerometer_Z)[2]), wrist_AccZ_fd,
  seq(0, 4, length=4), pace_fd)[,1]
  corfd_height_wristACCZ <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$wrist_Accelerometer_Z)[2], length=dim(input_data_fpca_30sec$wrist_Accelerometer_Z)[2]), wrist_AccZ_fd,
  seq(0, 4, length=4), height_fd)[,1]
  corfd_weight_wristACCZ <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$wrist_Accelerometer_Z)[2], length=dim(input_data_fpca_30sec$wrist_Accelerometer_Z)[2]), wrist_AccZ_fd,
  seq(0, 4, length=4), weight_fd)[,1]
  corfd_hip_wristACCZ <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$wrist_Accelerometer_Z)[2], length=dim(input_data_fpca_30sec$wrist_Accelerometer_Z)[2]), wrist_AccZ_fd,
  seq(0, 4, length=4), hip_fd)[,1]
  corfd_waist_wristACCZ <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$wrist_Accelerometer_Z)[2], length=dim(input_data_fpca_30sec$wrist_Accelerometer_Z)[2]), wrist_AccZ_fd,
  seq(0, 4, length=4), waist_fd)[,1]
  PC_no = seq(1,dim(input_data_fpca_30sec$wrist_Accelerometer_Z)[2],1)
  wristACCZ_corfd_summary = data.frame(PC_no,corfd_load_wristACCZ, corfd_pace_wristACCZ, corfd_height_wristACCZ, corfd_weight_wristACCZ, corfd_hip_wristACCZ, corfd_waist_wristACCZ )
  wristACCZ_corfd_summary_long = pivot_longer(wristACCZ_corfd_summary, cols = contains('wristACCZ') ,names_to = 'factors', values_to = 'corr_value')
  
  wristACCZ_corfd_plot = ggplot(data=wristACCZ_corfd_summary_long, aes(x= PC_no, y= corr_value))+
  geom_line() +
  geom_point()+ 
    facet_wrap(~factors, ncol =3)+
  ggtitle('Correlation coefficients of different factors and PCs of wrist ACC on Z direction')+ 
  ggplot2::labs(
    x = 'PCs', y = 'Correlation value') +
  ggplot2::theme(legend.position = 'bottom') +
  ggplot2::theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  )+
  ggplot2::theme(
     axis.title = ggplot2::element_text(size = 14),  # Increase axis label size
     axis.text = ggplot2::element_text(size = 14) # Increase axis text size
  ) 
  
  print(wristACCZ_corfd_plot)

##for Upperarm_acc_x
  
  corfd_load_upperarmACCX <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$upperarm_Accelerometer_X)[2], length=dim(input_data_fpca_30sec$upperarm_Accelerometer_X)[2]), upperarm_AccX_fd,
  seq(0, 4, length=4), load_fd)[,1]
  corfd_pace_upperarmACCX <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$upperarm_Accelerometer_X)[2], length=dim(input_data_fpca_30sec$upperarm_Accelerometer_X)[2]), upperarm_AccX_fd,
  seq(0, 4, length=4), pace_fd)[,1]
  corfd_height_upperarmACCX <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$upperarm_Accelerometer_X)[2], length=dim(input_data_fpca_30sec$upperarm_Accelerometer_X)[2]), upperarm_AccX_fd,
  seq(0, 4, length=4), height_fd)[,1]
  corfd_weight_upperarmACCX <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$upperarm_Accelerometer_X)[2], length=dim(input_data_fpca_30sec$upperarm_Accelerometer_X)[2]), upperarm_AccX_fd,
  seq(0, 4, length=4), weight_fd)[,1]
  corfd_hip_upperarmACCX <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$upperarm_Accelerometer_X)[2], length=dim(input_data_fpca_30sec$upperarm_Accelerometer_X)[2]), upperarm_AccX_fd,
  seq(0, 4, length=4), hip_fd)[,1]
  corfd_waist_upperarmACCX <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$upperarm_Accelerometer_X)[2], length=dim(input_data_fpca_30sec$upperarm_Accelerometer_X)[2]), upperarm_AccX_fd,
  seq(0, 4, length=4), waist_fd)[,1]
  PC_no = seq(1,dim(input_data_fpca_30sec$upperarm_Accelerometer_X)[2],1)
  upperarmACCX_corfd_summary = data.frame(PC_no,corfd_load_upperarmACCX, corfd_pace_upperarmACCX, corfd_height_upperarmACCX, corfd_weight_upperarmACCX, corfd_hip_upperarmACCX, corfd_waist_upperarmACCX )
  upperarmACCX_corfd_summary_long = pivot_longer(upperarmACCX_corfd_summary, cols = contains('upperarmACCX') ,names_to = 'factors', values_to = 'corr_value')
  
  upperarmACCX_corfd_plot = ggplot(data=upperarmACCX_corfd_summary_long, aes(x= PC_no, y= corr_value))+
  geom_line() +
  geom_point()+ 
    facet_wrap(~factors, ncol =3)+
  ggtitle('Correlation coefficients of different factors and PCs of upperarm ACC on X direction')+ 
  ggplot2::labs(
    x = 'PCs', y = 'Correlation value') +
  ggplot2::theme(legend.position = 'bottom') +
  ggplot2::theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  )+
  ggplot2::theme(
     axis.title = ggplot2::element_text(size = 14),  # Increase axis label size
     axis.text = ggplot2::element_text(size = 14) # Increase axis text size
  )
  
  print(upperarmACCX_corfd_plot)

  ##for Upperarm_acc_y
  
  corfd_load_upperarmACCY <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$upperarm_Accelerometer_Y)[2], length=dim(input_data_fpca_30sec$upperarm_Accelerometer_Y)[2]), upperarm_AccY_fd,
  seq(0, 4, length=4), load_fd)[,1]
  corfd_pace_upperarmACCY <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$upperarm_Accelerometer_Y)[2], length=dim(input_data_fpca_30sec$upperarm_Accelerometer_Y)[2]), upperarm_AccY_fd,
  seq(0, 4, length=4), pace_fd)[,1]
  corfd_height_upperarmACCY <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$upperarm_Accelerometer_Y)[2], length=dim(input_data_fpca_30sec$upperarm_Accelerometer_Y)[2]), upperarm_AccY_fd,
  seq(0, 4, length=4), height_fd)[,1]
  corfd_weight_upperarmACCY <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$upperarm_Accelerometer_Y)[2], length=dim(input_data_fpca_30sec$upperarm_Accelerometer_Y)[2]), upperarm_AccY_fd,
  seq(0, 4, length=4), weight_fd)[,1]
  corfd_hip_upperarmACCY <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$upperarm_Accelerometer_Y)[2], length=dim(input_data_fpca_30sec$upperarm_Accelerometer_Y)[2]), upperarm_AccY_fd,
  seq(0, 4, length=4), hip_fd)[,1]
  corfd_waist_upperarmACCY <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$upperarm_Accelerometer_Y)[2], length=dim(input_data_fpca_30sec$upperarm_Accelerometer_Y)[2]), upperarm_AccY_fd,
  seq(0, 4, length=4), waist_fd)[,1]
  PC_no = seq(1,dim(input_data_fpca_30sec$upperarm_Accelerometer_Y)[2],1)
  upperarmACCY_corfd_summary = data.frame(PC_no,corfd_load_upperarmACCY, corfd_pace_upperarmACCY, corfd_height_upperarmACCY, corfd_weight_upperarmACCY, corfd_hip_upperarmACCY, corfd_waist_upperarmACCY )
  upperarmACCY_corfd_summary_long = pivot_longer(upperarmACCY_corfd_summary, cols = contains('upperarmACCY') ,names_to = 'factors', values_to = 'corr_value')
  
  upperarmACCY_corfd_plot = ggplot(data=upperarmACCY_corfd_summary_long, aes(x= PC_no, y= corr_value))+
  geom_line() +
  geom_point()+ 
    facet_wrap(~factors, ncol =3)+
  ggtitle('Correlation coefficients of different factors and PCs of upperarm ACC on Y direction')+ 
  ggplot2::labs(
    x = 'PCs', y = 'Correlation value') +
  ggplot2::theme(legend.position = 'bottom') +
  ggplot2::theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  )+
  ggplot2::theme(
     axis.title = ggplot2::element_text(size = 14),  # Increase axis label size
     axis.text = ggplot2::element_text(size = 14) # Increase axis text size
  )
  
  print(upperarmACCY_corfd_plot)

  ##for Upperarm_acc_z
  
  corfd_load_upperarmACCZ <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$upperarm_Accelerometer_Z)[2], length=dim(input_data_fpca_30sec$upperarm_Accelerometer_Z)[2]), upperarm_AccZ_fd,
  seq(0, 4, length=4), load_fd)[,1]
  corfd_pace_upperarmACCZ <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$upperarm_Accelerometer_Z)[2], length=dim(input_data_fpca_30sec$upperarm_Accelerometer_Z)[2]), upperarm_AccZ_fd,
  seq(0, 4, length=4), pace_fd)[,1]
  corfd_height_upperarmACCZ <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$upperarm_Accelerometer_Z)[2], length=dim(input_data_fpca_30sec$upperarm_Accelerometer_Z)[2]), upperarm_AccZ_fd,
  seq(0, 4, length=4), height_fd)[,1]
  corfd_weight_upperarmACCZ <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$upperarm_Accelerometer_Z)[2], length=dim(input_data_fpca_30sec$upperarm_Accelerometer_Z)[2]), upperarm_AccZ_fd,
  seq(0, 4, length=4), weight_fd)[,1]
  corfd_hip_upperarmACCZ <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$upperarm_Accelerometer_Z)[2], length=dim(input_data_fpca_30sec$upperarm_Accelerometer_Z)[2]), upperarm_AccZ_fd,
  seq(0, 4, length=4), hip_fd)[,1]
  corfd_waist_upperarmACCZ <- fda::cor.fd(seq(0, dim(input_data_fpca_30sec$upperarm_Accelerometer_Z)[2], length=dim(input_data_fpca_30sec$upperarm_Accelerometer_Z)[2]), upperarm_AccZ_fd,
  seq(0, 4, length=4), waist_fd)[,1]
  PC_no = seq(1,dim(input_data_fpca_30sec$upperarm_Accelerometer_Z)[2],1)
  upperarmACCZ_corfd_summary = data.frame(PC_no,corfd_load_upperarmACCZ, corfd_pace_upperarmACCZ, corfd_height_upperarmACCZ, corfd_weight_upperarmACCZ, corfd_hip_upperarmACCZ, corfd_waist_upperarmACCZ )
  upperarmACCZ_corfd_summary_long = pivot_longer(upperarmACCZ_corfd_summary, cols = contains('upperarmACCZ') ,names_to = 'factors', values_to = 'corr_value')
  
  upperarmACCZ_corfd_plot = ggplot(data=upperarmACCZ_corfd_summary_long, aes(x= PC_no, y= corr_value))+
  geom_line() +
  geom_point()+ 
    facet_wrap(~factors, ncol =3)+
  ggtitle('Correlation coefficients of different factors and PCs of upperarm ACC on Z direction')+ 
  ggplot2::labs(
    x = 'PCs', y = 'Correlation value') +
  ggplot2::theme(legend.position = 'bottom') +
  ggplot2::theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  )+
  ggplot2::theme(
     axis.title = ggplot2::element_text(size = 14),  # Increase axis label size
     axis.text = ggplot2::element_text(size = 14) # Increase axis text size
  )  

print(upperarmACCZ_corfd_plot)

9 Showing all the plots of correlation coefficients of Acc signal vs. scalar inputs

trunkACCX_corfd_plot

trunkACCY_corfd_plot

trunkACCZ_corfd_plot

wristACCX_corfd_plot

wristACCY_corfd_plot

wristACCZ_corfd_plot

upperarmACCX_corfd_plot

upperarmACCY_corfd_plot

upperarmACCZ_corfd_plot

10 Plotting FPCs for Upperarm on X axis

df_ACCX_upperarm_2.5_15 = upperarm_total_df_pcs_30sec |> 
  dplyr::select(contains("Accelerometer_X"))|>
  dplyr::mutate(sub_con = rownames(upperarm_total_df_pcs_30sec)) |>
  pivot_longer(cols = contains('upperarm') ,names_to = 'PCs', values_to = 'FPCs_Values') |>
  dplyr::mutate(PCs = str_replace(PCs, ".*upperarm_Accelerometer_X_pc", "")) |>
  dplyr::filter(str_detect(sub_con,'2.5-15'))
df_ACCX_upperarm_2.5_15$PCs = as.integer(df_ACCX_upperarm_2.5_15$PCs)

upperarmACCX_FPCs_plot = ggplot(data=df_ACCX_upperarm_2.5_15, aes(x= PCs, y= FPCs_Values, group= sub_con))+
  geom_line() +
  geom_point()+ 
  facet_wrap(~sub_con, ncol =3, scales = "free")+
  ggtitle('FPCs of upperarm ACC on X direction for condition 2.5 kg load and 15 bpm pace')+ 
  #scale_color_grey()+
  ggplot2::labs(
    x = 'PCs', y = 'FPCs Scores') +
  ggplot2::theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  )+
  scale_x_continuous(breaks = seq(0, 12, by = 2))+
  ggplot2::theme(
     axis.title = ggplot2::element_text(size = 14),  # Increase axis label size
     axis.text = ggplot2::element_text(size = 14) # Increase axis text size
  )

print(upperarmACCX_FPCs_plot)

11 Residual Analysis for Functional Models

We implemented a goodness of fit check for our functional regression models via residual processes. Similar to the residual process of classical regression models, in the case of functional regression, we also seek to determine any dependencies between the residuals and predicted values that may indicate a lack of fit.

The code chunks below are to extract the fitted and residual values from our models of interest on which we then performed functional principal component analysis (FPCA) to obtain the FPCA scores for plotting.

## Remove unnecessary models
match <- c("_mag", "_acc2_", "_gyr2_", "_accgyr[x-z]_")
rm(list = c(unique(grep(paste(match,collapse="|"), ls(), value=TRUE))))

## Rename our models of interest
Task <- model_1
Task_individuals <- model_2
Individuals <- model_3

Stat_Mean_A_B <- model_mean
Stat_Rms_A_B <- model_rms
Stat_Med_A_B <- model_med
Stat_Perc75_A_B <- model_Perc75
Stat_Perc25_A_B <- model_Perc25
Stat_Min_A_B <- model_min
Stat_Max_A_B <- model_max
Stat_Mad_A_B <- model_mad
Stat_Std_A_B <- model_std


Med_A_B <- model_acc_sensors
Med_AG_W <- model_wrist_sensors
Med_AG_T <- model_trunk_sensors
Med_AG_U <- model_upperarm_sensors
Med_G_B <- model_gyr_sensors

FPC_G_B <- model_gyr_pcs_30sec
FPC_A_B <- model_acc_pcs_30sec
FPC_AG_W <- model_wrist_pcs_30sec
FPC_AG_T <- model_trunk_pcs_30sec
FPC_AG_U <- model_upperarm_pcs_30sec

rm(list = c(grep("model_", ls(), value = T))) 

all_models = list("Task and Individual Factors" = c("Task", "Individuals", "Task_individuals"),
                  "Statisitcal Feature Models" = c(grep("^Stat_", ls(), value = T)), 
                  "Median Feature Models" = c(grep("^Med_", ls(), value = T)), 
                  "FPC Models" = c(grep("^FPC_", ls(), value = T)))

# Custom function
fpca = function(x, name, threshold = 0.90) {
  ## prepare input
  x_matrix = matrix(x, ncol = 10, byrow = T) |> t()
  
  ## transform to a fd object
  argvals = seq(0,45,len=10)
  nbasis = 10
  basisobj = create.bspline.basis(c(0,45),nbasis)
  fdobj = smooth.basis(argvals, y=x_matrix, fdParobj=fdPar(basisobj))
  
  ## perform FPCA
  fpc = pca.fd(fdobj$fd, nharm=10)
  summed_variance = cumsum(fpc$varprop)
  num_fpc = which(summed_variance > threshold) |> min() ## determine the number of FPCs based on the threshold for CPV
  
  ## append all FPC scores to a df
  df = tibble(.rows = length(x)/10)
  cols = c()
  
  for (i in 1:num_fpc) {
  df = dplyr::bind_cols(df, fpc$scores[,i])
  cols = append(cols, paste0(name, "_fpc_", i))
  }
  
  colnames(df) = cols

  return(df)
}

11.1 Residual Plots

We constructed functional residual plots that comprise all pairwise plots of the residual FPC scores versus the fitted FPC scores to see if there are any discernible patterns or trends, similar to traditional residual plots.

Task and Individual Factors

Task

Individuals

Task_individuals

Statisitcal Feature Models

Stat_Mad_A_B

Stat_Max_A_B

Stat_Mean_A_B

Stat_Med_A_B

Stat_Min_A_B

Stat_Perc25_A_B

Stat_Perc75_A_B

Stat_Rms_A_B

Stat_Std_A_B

Median Feature Models

Med_A_B

Med_AG_T

Med_AG_U

Med_AG_W

Med_G_B

FPC Models

FPC_A_B

FPC_AG_T

FPC_AG_U

FPC_AG_W

FPC_G_B